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Abstract 

A  systematic  iterative  procedure  is  described  which 
yields  finite  difference  solutions  to  steady  free  boundary 
problems  of  plasma  physics  and  fluid  dynamics.   By  application 
of  the  method  of  steepest  descent  to  minimize  certain  energy 
functionals  connected  with  conformal  mapping  and  free  boundary 
problems,  a  system  of  artificially  time-dependent  equations 
is  produced  which  converges  to  a  steady  state  as  the  time 
parameter  goes  to  infinity.   The  solution  is  obtained  in  a 
parametric  form  in  which  the  physical  plane  is  the  conformal 
image  of  a  fixed  rectangle.   The  use  of  a  fixed  rectangle  as 
the  parameter  domain  enables  us  to  treat  effectively  the 
difficulties  engendered  by  the  singular  behavior  of  the  stream 
function  at  the  point  of  separation  between  the  fixed  and  free 
boundaries.   As  a  result  the  scheme  can  be  made  convergent  to 
second  order  in  the  mesh  width  h. 

The  method  has  been  applied  successfully  to  the  vena 
contracta  and  Riabouchinsky  models  of  classical  hydrodynamics 
for  both  plane  and  axially  symmetric  geometries.   An  accuracy 
of  four  significant  figures  in  the  calculation  of  quantities 
of  physical  interest  can  be  obtained  without  difficulty.   A 
plane  equilibrium  model  of  plasma  physics  is  also  computed. 
The  scheme  can  be  applied  to  a  variety  of  problems  with  two 
independent  variables. 


V 


1.  INTRODUCTION 

Free  boundary  problems  are  among  the  most  difficult 
to  treat  either  analytically  or  numerically.   Examples  of 
such  problems  are  found  in  the  occurrence  of  jets  or  cavities 
in  classical  hydrodynamics,  where  the  location  of  a  portion 
r  of  the  boundary  is  unknown.   These  problems  can  be  formula- 
ted in  terms  of  a  stream  function  ^  which  is  the  solution  of 
a  Dirichlet  problem  in  the  flow  region  of  the  physical  plane. 
The  additional  free  boundary  condition 

q  =  constant 

is  specified  on  the  free  boundary  T,    where  q  is  the  speed  of 
the  fluid.   This  condition,  which  is  motivated  physically  by 
Bernouilli's  Law,  compensates  for  the  fact  that  r  is  not 
known.   In  the  case  of  plane  flow,  where  function  theory  can 
be  applied,  the  solution  in  certain  cases  can  be  found 
analytically  by  the  hodograph  method.   For  three  dimensional 
flows,  for  example  those  with  axial  symmetry,  solutions  to 
familiar  models  are  not  known,  although  existence  and  unique- 
ness theorems  have  been  established.   Other  instances  of  free 
boundary  problems  which  are  of  great  interest  arise  in  plasma 
physics.   The  equilibrium  state  of  a  perfectly  conducting 
fluid  contained  by  a  magnetic  field  is  an  example,  where  the 
free  boundary  is  determined  by  the  balance  between  hydro-  . 
dynamic  and  magnetic  pressure. 


The  subject  of  this  paper  is  the  formulation  and 
implementation  of  a  finite  difference  method  for  treating 
these  problems  which  is  both  stable  and  consistent.   There 
are  two  difficulties  which  must  be  overcome  in  such  an 
approach: 

1)  It  can  be  shown  quite  generally  that  the  stream 

function  possesses  singular  behavior  at  the  separation  point 

z  between  the  free  boundary  r  and  the  fixed  boundary  Ft-, 
O  r 

of  the  flow  region.   More  precisely,  the  stream  function  has 

a  regular  expansion  there  in  terms  of  the  square  root 

1/2 
(z-z  )  '     ,    where  z  is  the  complex  variable  x +  iy  in  the 

physical  plane.   If  z  is  used  as  the  independent  variable, 

then  this  behavior  is  an  inherent  source  of  high  truncation 

error  in  the  discrete  solution. 

2)  Provision  must  be  made  for  the  determination  of  the 
free  boundary  by  some  systematic  iterative  procedure.   Again, 
if  the  physical  plane  is  used  as  the  domain  of  the  independent 
variable,  there  are  difficulties  associated  with  the  implemen- 
tation of  any  iterative  technique  caused  by  the  changing 
geometry  of  the  flow  region. 

The  method  we  shall  describe  is  designed  to  deal  with 
both  of  these  difficulties  and  can  be  applied  to  a  variety 
of  problems  involving  only  two  independent  variables.   We 
establish  a  parametric  formulation  by  introducing  a  fixed 
rectangle  R  in  an  auxiliary  w-plane.   We  seek  the  stream 
function  as  the  solution  to  an  appropriate  Dirichlet  problem 


in  R,  and  we  seek  the  flow  region  as  the  image  of  R  by  a 
conformal  mapping  z  =  z(w).   We  can  set  up  the  boundary 
correspondences  for  this  mapping  so  that  a  corner  of  the 
rectangle  is  mapped  into  the  separation  point.   When  this  is 
done,  the  difficulties  caused  by  singular  behavior  are  removed 
because  both  ip   and  z  become  regular  functions  of  w  at  the 
separation  point.   The  required  conformal  mapping  is  obtained 
by  solving  an  artificially  time-dependent  system  of  equations 
which  converge  to  a  steady  state.   This  system  of  equations 
is  derived  by  applying  the  method  of  steepest  descent  to 
minimize  certain  energy  functionals  connected  with  conformal 
mapping  and  free  boundary  problems.   Steady  state  is  achieved 
under  conditions  which  guarantee  that  the  mapping  is  conformal 
if  and  only  if  the  free  boundary  condition  is  satisfied  on  r. 
The  difficulty  of  locating  the  free  boundary  is  automatically 
eliminated  because  we  solve  all  equations  in  the  fixed 

rectangle  R. 

It  should  be  emphasized  that  since  we  have  uniformized 
the  singularity  at  the  separation  point  in  this  formulation, 
we  expect  the  over-all  finite  difference  method  to  converge 
to  the  continuous  solution  as  the  mesh  width  h  -^  0.   Further- 
more, since  all  equations  are  solved  in  a  rectangle,  it 
becomes  a  particularly  simple  task  to  write  finite  difference 
analogs  which  are  accurate  to  second  order  in  h.   Therefore 
the  convergence  is  O(h^).   One  of  the  objects  of  our  numerical 
calculations  is  to  verify  this  behavior  by  using  a  sufficiently 
large  number  of  mesh  points.   Our  use  of  conformal  mapping  is 


quite  essential  because  It  assures  that  the  equation  for  the 
stream  function,  when  written  In  a  conformally  invariant 
form,  makes  sense  In  the  rectangle.   This  conformal  Invariance 
Is  obvious  in  the  case  of  plane  flows,  but  it  exists  also, 
for  example,  in  the  case  of  axially  symmetric  flows  when  the 
governing  equation  is  written  in  the  form 

A^  -  (Vy.  V^)/y  =  0  . 

In  general,  the  modification  of  the  scheme  to  deal  with  any 
example  in  two  independent  variables  whose  differential  equa- 
tion is  conformally  invariant  is  a  relatively  easy  matter. 

The  central  issue  in  this  method  is  not  the  solution  of 
the  equation  for  f ,    as  it  is  in  the  direct  approach,  but  the 
establishment  of  the  appropriate  flow  region  as  the  conformal 
image  of  the  rectangle.   In  achieving  this  mapping  the  stream 
function  enters  only  in  the  determination  of  the  speed  along 
the  free  boundary.   For  example,  in  axially  symmetric  flows, 
the  speed  along  any  streamline  is  given  by 

q  =  |v^|/y  , 

and  a  conformal  mapping  is  sought  for  which  this  quantity  is 
constant  along  r.   We  include  calculations  of  axially  symmetric 
flows  in  our  sample  problems  to  demonstrate  the  applicability 
of  the  method  in  these  cases. 


In  Section  2  we  describe  the  particular  models  which 
were  studied  as  examples.   These  are  the  vena  contracta  and 
Rlabouchlnsky  models  of  classical  hydrodynamics  for  both  plane 
and  axlally  symmetric  geometries,  and  a  simplified  plane 
equilibrium  model  of  plasma  physics,  which  provided  a  pilot 
case  on  which  the  method  was  tested.   We  also  give  details 
relative  to  the  calculations  which  were  made  for  these  models, 
and  a  brief  review  of  previous  research  effort.   Section  3 
contains  the  derivation  of  the  method.   In  Section  4  details 
of  the  application  of  the  method  to  each  of  the  models  are 
described,  as  well  as  the  general  numerical  formulation.   The 
bulk  of  the  research  work  was  done  here,  for  there  are  several 
subtle  considerations  which  arise  in  connection  with  the 
problem  of  mapping  the  corner  of  R  into  the  separation  point. 
In  Section  5  we  give  the  results  of  the  computations.   Our 
fundamental  conclusion  is  that  the  method  can  be  applied 
successfully  to  the  models  we  studied,  although  In  some  cases 
convergence  is  a  delicate  matter.   Beyond  this  we  verify  that 
the  convergence  is  0(h  )  and  we  demonstrate  the  accuracy 
attainable  by  calculating  contraction  coefficients  and  drag 
coefficients  to  four  significant  figures.   Also,  a  few 
fragmentary  results  are  presented  pertinent  to  the  use  of  our 
method  as  a  technique  for  conformal  mapping  per  se.   Appendix 
B  comprises  the  F0RTRAN  IV  coded  program  for  both  the  plane 
and  axlally  symmetric  Rlabouchlnsky  models,  as  written  for 
the  CDC  6600  computer  at  New  York  University. 


2.  FREE  BOUNDARY  MODELS 

We  describe  the  three  mathematical  models  to  which  our 
numerical  method  was  applied.   In  each  case  a  brief  descrip- 
tion of  the  physical  problem  is  given,  the  equations  of 
steady  flow  are  formulated,  and  analytic  solutions,  where 
known,  are  given.   For  the  vena  contracta  model,  an  asymptotic 
analysis  of  the  behavior  of  the  jet  is  also  described,  to 
provide  a  theoretical  basis  for  certain  extrapolations  used 
in  calculating  contraction  coefficients.   We  close  with  a 
brief  discussion  of  some  of  the  previous  research  on  these 
problems  and  the  difficulties  inherent  in  a  finite  difference 
approach. 

2.1.  Cusped  Geometry  Model 

Our  purpose  in  studying  this  model  is  that  it  provides 
an  elementary  example  with  which  to  exhibit  the  numerical 
technique  to  be  described.   It  incorporates  some,  but  not  all, 
of  the  features  of  the  method.   The  cusped  geometry  model  is 
of  importance  in  plasma  physics  as  a  simplified  plane  mathe- 
matical model  of  the  equilibrium  state  of  a  perfectly  con- 
ducting fluid  contained  by  a  magnetic  field  [1]. 

The  fluid  is  located  in  a  doubly  periodic  array  of  cusped 
figures  of  equilibrium  with  boundaries  F,  whose  cross-sections 
are  shown  in  Figure  2.1.   Their  centers  lie  at  the  lattice 
points  (n+in')-^  (n  and  n'  odd).   A  magnetic  field  B  is  genera- 
ted by  currents  which  flow  in  alternate  directions  along 
circular  coils  of  radius  r  and  boundary  O  which  are  also 


periodically  arranged  at  the  points  (n+lnM-^  (n  and  n'  even), 
In  the  vacuum  region  D  between  r  and  O,  the  field  B  satisfies 
the  Maxwell  equations 

V  XB  =  0 

V  .  B  =  0 

with  the  associated  boundary  condition 

B.  V  =  0 

on  O  and  r.   Since  the  vector  B  lies  in  the  (x,y)-plane  we 
can  Introduce  a  scalar  magnetic  potential  ip   defined  by 

.  B  =  V  x(0,0,t) 

and  formulate  the  problem  for  f,    after  Introducing  an  appro- 
priate normalization,  as 

A^  =  0 
in  D, 

^  =  0 
on  r,  and 

f  =  ±1 

on  the  boundary  O  of  alternate  coils.  The  balance  between 
magnetic  forces  and  the  hydrodynamic  pressure  in  the  fluid 
gives  the  additional  free  boundary  condition 


2       2 

B  =  (V^)   =  constant 


which  holds  on  r. 

The  uniqueness  of  the  solution  and  the  stability  of 
these  cusped  figures  of  equilibrium  can  be  established  by 
showing  that  the  energy  of  the  system  is  a  minimum  relative 
to  all  configurations  with  the  same  total  mass  of  fluid  [1  ]. 

For  our  numerical  method,  it  is  convenient  to  formulate 
this  problem  in  a  different  way.   Considerations  of  symmetry 
permit  us  to  confine  our  attention  to  the  determination  of  tl/ 
in  the  fundamental  region  D  bounded  by  the  lines  x  =  ±£, 
Y  =  ±£   and  the  corresponding  parts  of  r,  as  shown  in 
Figure  2.1.   By  symmetry,  we  have 


f  =   0 


along  these  lines.   We  introduce  a  fixed  rectangle  R  in  the 
first  quadrant  of  an  auxiliary  w-plane,  where  w  =  u+iv.   The 
rectangle  is  bounded  by  the  coordinate  axes  and  the  lines 
u  =  a,  V  =  b.   We  seek  the  solution  parametrically  in  terms 
of  the  variable  w  in  R.   The  mapping  from  R  to  D   is  to  be  a 
conformal  transformation,  periodic  in  u,  in  which  the  boundary 
V  =  b  maps  onto  O,  the  boundary  v  =  0  maps  onto  the  outer 
boundary  of  D  ,  and  the  vertical  boundaries  u  =  0,  u  =  a  map 
onto  the  segment  y  =  0,  r  <_  x   <_  £.      With  these  boundary 
correspondences,  the  parametric  solution  for  ij/   is  given  by 


^  =  v/b 
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when  the  boundary  conditions  on  bD     are  transferred  to  the 
boundary  of  R.   Apart  from  the  factor  b,  we  observe  that  R 
Is  the  domain  of  a  complex  potential  of  which  f   is  the 
imaginary  part.   Our  original  problem  is  therefore  inverted, 
and  becomes  the  determination  of  the  unique  conformal  mapping 
of  R  onto  D  ,  with  the  indicated  boundary  correspondences, 
for  which  the  boundary  r  is  located  so  that  the  free  boundary 
condition  holds  along  it.   Relative  to  the  w-plane,  the  free 
boundary  condition  takes  the  form 


|V^^|   =  tp^/ix^  +  Y^)    =   constant  . 


The  technique  for  achieving  this  conformal  mapping  is  the 
basis  of  the  numerical  method  and  is  described  in  Section  3» 

2.2.   Riabouchinsky  Model 

This  hydrodynamlc  model  is  a  representation  of  steady 
potential  flow  past  an  obstacle  C  in  which  cavitation  occurs 
[  2  ].   The  cavity  volume  is  kept  finite  by  the  introduction 
of  a  fictitious  image  obstacle  C*  located  downstream  from  C, 
We  study  this  flow  in  both  plane  and  axially  symmetric 
geometries  for  the  configuration  shown  in  Figure  2.2. 

A  uniform  stream  directed  along  the  positive  x-axis  and 
of  speed  U  impinges  on  the  obstacle  C.   A  cavity  of  water 
vapor  is  formed  behind  C,  bounded  by  the  free  streamlines  r 
and  r*  which  detach  from  C  and  join  the  image  obstacle  C*  at 
the  points  indicated.   The  coordinates  are  placed  so  that 


symmetry  is  achieved  in  both  the  x-  and  y-axes.   In  the  plane 
case  these  coordinates  lie  In  a  cross-section  of  the  flow  and 
the  obstacle  Is  an  Infinite  flat  plate  whose  cross-section  C 
has  length  2y  .   In  the  axlally  symmetric  case.  Figure  2.2 
depicts  a  meridian  plane,  and  the  obstacle  Is  a  flat  circular 
plate  of  radius  y  .   The  flow  for  both  geometries  can  be 
described  In  terms  of  a  stream  function  t(/   =  i/{x,y),   defined  In 
the  region  D  exterior  to  the  cavity,  from  which  the  velocity 
components  of  a  fluid  particle  are  found  by  the  equations 


V  =  11/  /y 
X    V 


,m 


(1) 


V  =  -11/  /y 


m 


where  m  =  0  in  the  plane  case  and  m  =  1  in  the  axlally  symmetric 
case.   From  (1)  we  derive  the  relation  for  the  speed  q  along 
any  streamline  ^  =  constant,  as 


m    ,  /.m 


(2)  q  =  \^^^\/l     =  f^/y 


where  n  is  the  normal  to  the  streamline  In  the  usual  sense, 
and  V  denotes  the  gradient  in  the  z -plane,  z  =  x+iy.   The 
free  streamlines  r  and  F  are  unknowns  and  must  be  determined 
as  part  of  the  solution.   The  cavity  is  a  domain  of  constant 
vapor  pressure,  and  for  stability  of  the  flow  we  expect  a 
continuous  variation  of  pressure  across  the  free  streamlines. 
Therefore,  an  application  of  Bernouilli's  law  yields  the  free 
boundary  condition 


10 


q  =  q  =  constant 


along  r  and  r*.   We  define  the  cavitation  parameter  o   by  the 
formula 

The  significance  of  a  Is  that  specification  of  Its  value  fixes 
the  shape  of  the  cavity,  apart  from  a  multiplicative  factor 
which  scales  the  size  of  the  physical  domain  D. 

Using  symmetry,  we  consider  the  flow  in  the  subdomain 
D  of  D  in  the  first  quadrant  of  the  z-plane.   We  have  the 
following  problem  for  the  determination  of  ^: 

(3a)  A^  -  (m/y).^y  =  0   in  D^  , 

{3b)  f  =  0     on     r+C*+x-axls    , 

(5c)  ^x  =  ° 

on  the  y-axls  above  the  free  streamline,  by  symmetry, 

(3d)  q  =  ^o   on  r  , 


^o 


{3e)         t  ~  U/(l+m)-y-'-'^"^  as   |  z  |  ^  oo 


The  asymptotic  behavior  (3e)  is  consistent  with  the  assumption 
of  a  uniform  stream  at  oo  . 


11 


In  the  plane  model,  where  we  have  the  tool  of  function 
theory,  the  solution  can  be  found  analytically  by  the  hodo- 
graph  method  [2],  in  terms  of  the  complex  potential  ^  =  ({)+  If, 
as  the  conformal  mapping 


^o 


c  /C^-  -^—   ^/^ 


a2(^2-l)V2^(l^a2)V2/,.|  (-^^^ 


d^ 


-I 


l/d-m^)'/" 


where  the  quantity  a  is  defined  by 


a  =  I  {%--a^)/%^ 


and  the  C-plane  has  been  normalized  so  that  C  =  1  at  the 
stagnation  point  B.   The  domains  D-  and  D  relative  to  {k)   are 
shown  in  Figure  2.3a  and  2.3b. 

For  the  axially  symmetric  case,  the  analytic  solution  is 
not  known.   Existence  and  uniqueness  proofs  have  been  estab- 
lished for  a  wide  variety  of  such  axially  symmetric  flows  by 
use  of  variational  principles  [3]. 

It  is  necessary,  for  our  numerical  technique,  to  express 
the  plane  solution  in  a  certain  parametric  form.   We  choose 
to  normalize  the  flow  so  that  q  =  1,  and  we  define  the 
parameter  k  by 


12 


k  =  (l+a2)-l/2 


We  Introduce  a  rectangle  R  In  a  w-plane,  w  =  u+iv,  oriented 
as  shown  In  Figure  2.3c,  and  observe  that  the  transformation 

(5)  C-k.sn{w,k) 

maps  R  onto  D^  [^  ].   The  ratio  k'/K  is  determined  uniquely  by 
k,  and  the  scaling  is  chosen  so  that  the  boundary  corre- 
spondences are  as  shown.   In  this  context,  the  parameter  k  is 
the  modulus  of  the  Jacobian  elliptic  function  sn.   Substitution 
of  (5)  into  the  mapping  (4),  and  use  of  the  normalization  we 
have  indicated,  gives 

(6)  z(w)  =  (l/k)[E(w,k)  -  k'^w +ik'dn(w,k)  +  i(E' -  k^K' )]  , 

where  E(w,k)  is  the  elliptic  Integral  of  the  second  kind, 

E'  =  E(K  ,k  ),  and  k'  is  the  complementary  modulus  defined  by 


1  '2      n     1  2 

k   =  1  -  k 


The  drag  coefficient  C^.,  the  cavitation  parameter  a,  and  the 
speed  at  infinity  U  can  be  expressed  in  terms  of  k  by  the 
formulas 

Cp  =  2{l+a)(E'  -k^K')/(k'^+E'  -  k^K*) 
a  =  2kV(l-k') 
U  =  k/(l+k')   . 
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These  quantities,  and  equations  (5)  and  (6),  are  used  both  to 
provide  initial  conditions  for  the  numerical  calculation  and 
to  check  the  accuracy  of  the  finite  difference  solution  in  the 
plane  case.   We  note  that  as  the  length-to-width  ratio  r  of 
the  rectangle  increases,  the  parameter  k  — »•  1 .   The  model 
approaches  classical  Helraholtz  flow  [  2  ]  in  which  o  —*■  0   and 
the  cavity  behind  the  plate  becomes  infinite.   In  the  limiting 
case,  we  have 

C^   =  27r/(Tr+4)  =   .87980 
U  =  1  . 

In  our  numerical  calculations  we  shall  estimate  the  drag 
coefficient  for  both  plane  and  axially  symmetric  cases  as  a  -*  0. 
It  should  be  observed  that  letting  r  -♦■  00  is  a  poor  way  of 
approximating  Helmholtz  flow,  for  it  is  well  known  that  the  shape 
of  the  Helmholtz  free  boundary  approaches  a  parabola  at  large 
distances  from  the  plate. 

2.3.  Vena  Contracta  Model 

This  is  a  hydrodynamlc  model  of  steady  potential  flow  of 
an  incompressible  fluid  which  issues  from  an  orifice  in  a  plane 
wall  [2].   The  geometry  is  shown  in  Figure  2.4  for  both  plane 
and  axially  symmetric  geometries.   The  physical  coordinates  x 
and  y  define  a  cross-section  of  the  flow  in  the  plane  case  and 
the  orifice  is  an  infinite  slit  of  width  2y^.   In  the  axially 
symmetric  case,  these  coordinates  depict  a  meridian  plane  and 
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the  orifice  is  to  be  Interpreted  as  a  circular  hole  of  radius 

^o- 

A  Jet  of  fluid  issues  horizontally  from  the  orifice,  to 

the  left  of  which  lies  an  infinite  reservoir  of  fluid.   The 

flow  can  be  characterized  in  both  geometries  by  a  stream 

function  -p  =   ■^(x,y)  from  which  the  velocity  components  and  the 

speed  are  given  by  equations  (1)  and  (2).   The  x-axis  is  a  line 

of  symmetry  so  that  the  flow  region  D   is  bounded  by  the  x-axis, 

by  the  vertical  wall  y  ^  y^  and  by  the  free  streamline  r  which 

detaches  from  the  wall  at  the  separation  point  z   =  iy  and  is 

^         0       0 

asymptotic  to  the  line  y  =  y   as  x  -*-+oo.   Along  this  free 
boundary,  whose  location  must  be  found  as  part  of  the  solution, 
continuity  of  the  pressure  and  an  application  of  Bernouilli's 
law  gives  the  condition  that  the  fluid  speed  there  is  constant. 
Thus  we  have  the  following  problem  for  the  determination  of  ij/: 

(7a)  A^  -  (m/y).^y  =  0   in  D^  , 

where  m  =  0  in  the  plane  case  and  m  =  1  in  the  axially  symmetric 
case, 

(7b)  if/  =  0        on  X-axis  , 

(7c)  f  =   C 

on  y  >  y  and  on  r,  where  C  7^  0  is  a  constant,  and 


(7d)  q  =  q    on  r  , 

^o 
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where  Q^^  is  the  constant  speed.  Par  downstream  in  the  jet 
the  flow  becomes  uniform,  so  that  f   has  the  asymptotic  form 


(7e)  t  ~  q^y-^'^'^/Cl+m)   as  x  ^  +00 


In  terms  of  the  speed  q  .   The  constant  c  In  (7c)  is  related 


to  the  other  parameters  by 


0  =  qo^Jc™/'!™) 


The  quantitative  measure  of  the  vena  contracta  (the  narrowing 

of  the  jet  as  X  — ►  +00  )  is  the  contraction  coefficient  C  defined 

m 

as  the  ratio  of  the  cross-sectional  area  of  the  jet  at  x  =  +00 
to  its  area  at  the  orifice.   For  the  geometries  we  study  this 
definition  takes  the  form 


c  =  (y  /y  )-^^'" 


In  the  plane  case  the  flow  can  be  found  analytically  by 
the  hodograph  method  in  terms  of  the  complex  potential  C  =  <t)  +  i^- 
Introducing  the  normalization  C  =  1,  we  obtain  the  solution  as 
the  conformal  map  of  the  strip  0  j<  Im  C  ^  1  onto  D  given  by 


(8)       z(C)  =  -(2C^yyTr) 


.-C7r/2^(,^,-7rC)l/2 


+  log  (e^    [ (1+e   ^)  ^   -  IJ)   , 
with  boundary  correspondences  arranged  so  that 
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z(i)  =  y^ 

z(+00  )  =  +00 

z(-oo  )  =  -oo  . 

By  substitution  of  ^  =  i  into  (8),  the  contraction  coefficient 
is  found  to  be 

Cq  =  7r/(Tr+2)  . 

The  analytic  solution  for  the  axially  symmetric  vena 
contracta  is  not  known. 

For  our  numerical  approach  it  is  necessary  to  consider 
a  modified  vena  contracta  model,  in  order  to  conveniently 
obtain  a  parametric  formulation  in  which  the  physical  plane  is 
a  conformal  mapping  from  a  fixed  rectangle  R  in  an  auxiliary 
w-plane.   To  achieve  this  modification,  we  truncate  the 
infinite  jet  at  a  distance  &   downstream  from  the  fixed  wall, 
and  impose  the  boundary  condition 

(7e')  V'x/y"'  =  0 

along  X  =  £,  which  means  physically  that  the  velocity  vector 
there  is  horizontal.   All  other  conditions  (7a-7d)  remain  the 
same.   This  modification  is  tantamount  to  the  assumption  of  an 
image  reservoir  of  fluid,  at  distance  2i   downstream,  back  into 
which  the  fluid  flows.   In  this  respect  the  truncated  vena 
contracta  model  is  quite  similar  to  the  Riabouchlnsky  model. 
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We  conjecture  that  the  flow  corresponding  to  the  trunca- 
ted jet  problem  converges  to  the  infinite  jet  flow  in  any- 
fixed  bounded  subdomain  of  D  to  the  left  of  x  =  ^  as  ^  -♦  oo  . 
In  particular,  we  define  the  truncated  jet  contraction 

coefficient  C  {£)    to  be  the  ratio  of  the  cross-sectional  area 
m 

of  the  jet  at  X  =  i  to  its  area  at  the  orifice,  and  we  assert 
that  the  asymptotic  behavior  of  C  (i)  is  given  by 

-A  & 

(9)  C  (^)  ~  C  +A  e  "^    as  ^  -»•  oo 

-^  m       mm 


in  both  plane  (m  =  0)  and  axially  symmetric  (m  =  1)  geometries, 
for  appropriate  values  of  the  constants  A  and  A  .   Let 
z  =  z(w)  be  the  analytic  function  we  seek  which  maps  the 
rectangle  R  in  the  w-plane  onto  the  domain  D  .   Assume  this 
mapping  is  arranged  so  that  the  boundary  correspondences  are  as 
shown  in  Figure  2.5,  and  is  normalized  so  that  the  jet  height 
at  X  =  ^  is  1.   Then  for  large  I   the  transformation  has  the 
form 


z(w)  =  a^  +  ttgW  +z^(w)  , 


where  a^  and  a^   are  real  constants  and  z.(w)  -*  0  as  Re  (w)  — »■  a; 
that  is,  the  mapping  is  approximately  affine  near  the  edge  BC. 
Prom  this  observation  we  conclude  that  I   satisfies  the 
approximate  formula 


I   =  a  +  r  , 
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where  r  Is  the  length-width  ratio  of  R.   Therefore  the 
asymptotic  formula  (9)  is  equivalent  to 


-A  r 
(10)        ^m^^^  '^  ^m'*"^m®         as  r  ^  oo 


for  some  constants  A_  and  A_. 

m      m 

In  what  follows  we  aim  first  to  directly  verify  (10)  in 
the  plane  case  and  second  to  obtain  estimates  for  the  values 
of  A  for  both  the  plane  and  axially  symmetric  models.   Our 
purpose  in  doing  this  is  to  provide  a  check  on  the  estimation 
of  C  by  our  numerical  method.   In  the  plane  case,  the  trunca- 
ted model  can  be  solved  completely  by  the  hodograph  method. 
Prom  this  solution  the  contraction  coefficient  can  be  calculated 
directly  as 

C^(r)  =  (k'+E')/(l  +  E') 

where  k'  is  the  complementary  modulus  of  the  elliptic  functions 
and  integrals  associated  with  the  rectangle  R,  and  E  is  the 
complete  elliptic  integral  of  the  second  kind  with  respect 
to  k'.   This  integral  is  even  in  k'  and  has  the  expansion  [^  ] 

•   e'  =  7r/2  +  0(k'^) 

The  parameter  k'  is  related  to  r  by  the  formula  [  4  ]  • 

k  =  e   '^   +  0{e    ) 
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for  large  r.   Therefore  we  obtain 


C    (r)   ^  ir/{-JT+2)  +A^e   ^^^      ,        as   r -♦  oo 


which  verifies  (10)  and  gives 

^o  =  ^/2 

for  the  plane  case. 

For  the  axially  symmetric  case,  where  we  have  no  rigorous 
analysis,  we  can  show  that  the  asymptotic  behavior  of  the  free 
streamline  for  the  infinite  jet  model,  where  y   is  normalized 


to  be  1,  is  given  by 

(11)       y(x)  ~  1  +  A^e  "   ,   X  -.  CO 


-An  X 


with  an  estimate  on  the  value  of  A^  which  appears  in  the 
exponent.   We  make  the  further  conjecture  that  this  value  of 
A,  is  the  same  constant  which  appears  in  (10).  To  prove  (11), 
we  observe  that  by  (7e)  the  stream  function  for  the  infinite 
jet  model  has  the  form 

(12)  f  =  y^/2  +   ^3_(x,y) 


whe 


re  ^n  — »■  0  as  X  — ►  +oo  ,  and  where  we  have  adopted  the 


normalization  q  =1.   The  function  ^  satisfies 
^o 
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Tj/   =  1/2 

on  the  free  boundary  r,  and  therefore  also  the  mixed  condition 

(13)  f^/y  +   {fc/y){f  -   1/2)  =  1  , 

where  K   is  the  curvature  along  r.   Following  an  idea  due  to 
Garabedian  [  5  ],  we  observe  that  the  condition  that  the  flow 
is  irrotational. 


q^  +  Kq  =  0 


which  holds  along  any  streamline,  can  be  used  to  show  that  (13) 
is  satisfied  by  ip   along  an  approximation  r*  to  r  to  second 
order  in  the  normal  displacement  5n  between  r  and  r*.   We  use 
this  result  to  obtain  an  estimate  for  A,.   Let  the  line  y  =  1 
be  the  approximate  free  boundary  r*.   We  impose  the  boundary 
condition  (13)  on  r*,  and  use  equation  (12)  and  the  fact  that 
K    =  0  on  r*  to  obtain  the  condition 


(tl)y  =  0 


on  r  .   On  the  x-axis  we  have 


f  =  f^  =   0   , 


and  throughout  D  ,  ^  satisfies 
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^^1 "  ^hK^^  =  °  • 


Thus  we  have  an  eigenvalue  problem  for  the  determination  of 
fj,    which  is  valid  for  large  positive  x.   Using  separation  of 
variables  we  obtain  the  solutions 


f^  =   Cye"^''.  J^(Ay) 


The  smallest  value  of  A  which  satisfies  the  boundary  conditions 
is  the  smallest  solution  to 


Jq(A)  +AJJ(A)  =  -A^  J^(A)  =  0  , 


from  which  we  obtain 


A^  =  2.4 


Substituting  f^    into  (12)  and  using  the  condition  f  =   l/2  on  r 
gives  an  implicit  equation  for  the  free  boundary 


1/2  =  Y^/2   +   Ce-^^'y  J^CAy)  . 


It  is  easy  to  show  that  for  sufficiently  large  x  this  equation 
has  the  form  (11),  which  is  what  we  set  out  to  prove.   The  same 
analysis  can  be  carried  out  for  the  plane  model.   The  result  is 
-^     =   7r/2  which  agrees  with  the  value  obtained  by  the  rigorous 
treatment.   This  supports  our  conjecture  that  the  free  boundary 
of  the  infinite  jet  and  the  truncated  contraction  coefficient 
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have  quantitatively  the  same  asymptotic  behavior.   We  use  these 
results  in  Section  4.5.7.   We  observe  here  that  the  value  of 
the  exponent  A  is  substantially  larger  in  the  axially 
symmetric  case  compared  to  the  plane  case.   This  suggests  that 
the  convergence  to  the  infinite  jet  model  as  i  -+  oo  is  sub- 
stantially faster  in  the  axially  symmetric  case  than  in  the 
plane  ease. 

2.4.  The  Difficulties  of  Applying  a  Finite-Difference  Method 

For  any  of  the  models  introduced  there  are  two  problems 
associated  with  using  finite  differences  to  solve  directly  for 
^  as  a  function  of  the  physical  variables.   First,  the  deter- 
mination of  the  free  boundary  requires  an  Iterative  technique 
in  which  a  trial  free  boundary  is  selected  and  then  modified 
in  the  next  iterative  cycle  to  produce  a  better  guess.   Such 
schemes  have  been  applied,  with  the  corrections  made  essentially 
by  trial  and  error,  in  [ 6  ]  for  the  axially  symmetric  vena 
contracta,  and  in  [  7 ]  for  the  plane  and  axially  symmetric 
Riabouchinsky  models.   Corrections  of  a  trial  free  boundary 
can  be  made  automatically  with  the  technique  used  in  Section 
2.3  for  obtaining  the  asymptotic  behavior  of  the  vena  contracta 
Jet,  and  should  produce  a  second  order  rate  of  convergence. 
This  has  been  accomplished  recently  in  [8  ],  where  the  free 
boundary  is  approximated  by  choosing  the  multiplicative 
constants  in  a  linear  combination  of  appropriate  functions. 
The  second,  and  more  serious,  problem  is  that  the  separation 
point  z  is  a  singular  point  of  the  flow.   For  the  plane  models 
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the  hodograph  method  shows  that  the  stream  function  expansion 

1/2 
near  z  Is  regular  In  powers  of  (z-z^)  '    ,    and  the  same  type 

of  behavior  holds  in  the  axlally  symmetric  cases,  as  shown 

In  [  9 ] .   In  these  expansions  the  first  term  Is  zero  since  the 

speed  Is  finite  at  z  ,  but  In  any  case,  the  finite  difference 

equations  exhibit  unbounded  truncation  error  near  z  .   In  the 

papers  [  6  ]  and  [  7  ]  cited  above,  an  attempt  Is  made  to  deal 

with  this  difficulty  by  refinement  of  the  mesh  In  the  vicinity 

of  z  .   Such  procedures  localize  but  do  not  remove  the 
o 

difficulty.   The  results  of  these  calculations  are  not  con- 
vincing.  In  the  vena  contracta  model  the  free  boundary  makes 
an  angle  of  almost  30  with  the  vertical,  whereas  It  Is 
supposed  to  be  0°.   Also,  comparison  between  computed  and  exact 
Rlabouchinsky  plane  flows  shows  a  wide  discrepancy  In  the  free 
boundary  location.   In  [  8 ]  the  effect  of  the  singularity  Is 
reduced  by  Incorporating  Into  the  linear  combination  functions 
which  describe  the  singular  behavior.   This  makes  possible  a 
better  description  of  the  free  boundary  near  the  separation 
point. 

The  method  presented  in  this  paper  removes  both  of  these 
difficulties  in  its  application  to  the  hydrodynamlc  models.  We 
establish  a  systematic  iterative  procedure  for  the  determination 
of  r  which  is  made  independent  of  the  changing  geometry.  We 
deal  with  the  singular  behavior  by  formulating  the  problem 
parametrically.   The  physical  plane  is  sought  as  a  conformal 
transformation  from  a  fixed  rectangle,  one  of  whose  corners  is 
made  to  map  onto  the  separation  point.   We  describe  this  method 
in  detail  in  Section  3« 

2h 


5.  THE  METHOD  OP  STEEPEST  DESCENT 

3.1.  General  Description 

The  numerical  method  alters  the  formulation  of  a  free 
boundary  problem  in  the  following  way.   We  Introduce  a  fixed 
rectangle  R  in  an  auxiliary  w-plane  and  we  seek  the  solution 
in  the  parametric  form 

z  =  z(w)  . 

The  function  z(w)  is  to  be  found  as  the  analytic  function 
which  maps  R  onto  the  flow  region.   The  free  boundary  r  is  to 
be  determined  along  with  this  mapping.   The  function  ^'solves 
a  Dirichlet  problem  in  R,  provided  that  its  governing  equation 
is  conformally  invariant. 

The  possibility  of  implementing  this  approach  rests  on 
devising  a  technique  for  obtaining  the  conformal  map.   We 
couple  two  distinct  variational  principles,  one  for  free 
boundary  problems,  and  the  other  for  the  Plateau  problem,  which, 
by  the  application  of  the  method  of  steepest  descent,  produce 
an  artificially  time-dependent  system  of  equations  whose  steady 
state  solutions  yield  the  desired  mapping. 

3.2.  Variational  Principle  —  Conformal  Mapping 

We  cite  a  result  connected  with  the  minimum  formulation 
of  the  Plateau  problem  [lO].   Let  (^  be  the  class  of  differen- 
tiable  functions  z(w)  =  x(u, v) + ly(u, v)  which  map  the  domain  R 
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In  the  w-plane  onto  some  fixed  domain  D  in  the  z -plane,  and 
which  map  the  boundary  SR  onto  the  boundary  SD  in  such  a  way 
that  three  specified  points  on  SR  go  into  three  specified 
points  on  SD.   The  domain  R  is  a  rectangle  in  the  first 
quadrant  of  the  w-plane  with  horizontal  and  vertical  sides. 
Then  the  conformal  mapping  between  R  and  D  is  that  function 
in  G  which  minimizes  the  Dirichlet-Douglas  integral 

(1)  J  =  ^ff   [{Vx)2+  (Vy)2]dudv  . 

R 

The  three  point  condition  in  our  hypothesis  serves  to  make  the 
mapping  unique.   The  method  of  steepest  descent  can  be  applied 
to  J  in  order  to  find  this  conformal  mapping.   We  introduce  t 
as  an  artificial  time  parameter  and  regard  z  as  a  function  of 
t  as  well  as  w.   Then  J  becomes  a  function  of  t  which  we  want 
to  decay  as  rapidly  as  possible  as  t  -►  oo .   To  achieve  this, 
we  differentiate  J  with  respect  to  t  to  obtain 


=   -  J]    [Ax  x^  +  Ay  y^ldudv  -  ^0  (x^x^  +  y^y^)d£ 


dJ 
d^ 

R  SR 


where  v  is  the  inner  normal  and  s  is  the  arc  length  on  ^R. 
Paths  of  steepest  descent  are  those  which  assure  that  each  of 
the  integrands  on  the  right  is  positive.   Thus  we  are  led  to 
specify 

X,  =  aAx 
y^  =  aAy 
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in  R,  where  the  positive  constant  a  is  arbitrary.   This  factor 
controls  the  rate  of  convergence  as  t  — »■  oo  and  has  the  usual 
interpretation  in  the  steepest  descent  idea,  namely  that  the 
method  tells  the  "direction"  in  which  to  move,  but  not  how  far 
to  move.   On  the  boundary,  the  same  considerations  lead  to  the 
equation 

x^+iy^  -  P[(x^x^  +  y^y3)/(x^+y^)](x3+iy^)  . 

This  boundary  condition  is  consistent  with  the  requirement  that 
the  boundary  points  of  R  go  into  boundary  points  of  D.   It  has 
a  simple  geometric  interpretation.   The  vector  (x  ,y  )  is 
always  tangential  to  SD,  but  the  vector  (x  ,y  )  will  be  normal 
to  SD  only  when  steady  state  (x  x  +y  y  =  0)  is  reached.   A 

V  S    V  s 

boundary  point  is  therefore  shifted  tangentially  along  SD  at 
time  't  by  an  amount  proportional  to  the  tangential  component 
of  (x  ,y  ).   The  factor  p  is  a  positive  constant. 

As  t  -»■  00  a  steady  state  is  reached  in  which  we  obtain 

(2)  Ax  =  0  ,    Ay  =  0 
in  R,  and 

(3)  ^3+^3  =  ° 

on  Sr.   Since  the  rectangle  R  has  horizontal  and  vertical  sides 
the  s  and  v  derivatives  are  equivalent  to  u  and  v  derivatives, 
and  so  the  boundary  condition  (3)  implies 
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(k)  xx+vy=0 


on  SR.   The  function  x  x  +y  y  is  harmonic  in  R,  because  x 
and  y  are  themselves  harmonic  In  the  steady  state,  and  so  (4) 
actually  holds  throughout  R.   Moreover,  using  a  familiar 
theorem  from  function  theory,  we  obtain  the  result  that  the 
harmonic  conjugate  of  this  function  Is  Identically  constant, 
that  Is, 

t^y  2  ,  2     2    2   „ 

(5)  x^  +  y^  -  x^-y^  ^H 

In  R.   In  these  terms,  the  assertion  of  the  minimum  principle 
Is  that  H  =  0,  from  which  we  see  that  (4)  and  (5")  together  are 
equivalent  to  the  Cauchy-Rlemann  conditions.   The  method  of 
steepest  descent  therefore  provides  a  technique  for  obtaining 
conformal  mappings. 

Suppose  now  that  the  requirements  on  the  mapping  are 
over-specified  by  having  a  four-point  condition  instead  of  the 
three-point  condition  described  above.   If  for  our  models  we 
envision  fixed  guesses  for  the  free  boundary  locations,  then 
we  are  dealing  with  such  an  over-specified  problem.   If  the 
integral  (1)  possesses  a  minimum  at  all,  then  our  analysis  can 
be  carried  out  just  as  above  except  that  the  constant  H  in  (5) 
is  not  necessarily  equal  to  zero.   Thus  we  must  modify  the 
mapping  procedure  by  building  in  a  mechanism  for  moving  the 
free  boundary  r  around  until  in  the  steady  state  we  have  H  =  0. 
There  are,  of  course,  many  boundary  locations  for  which  this 
conformality  condition  is  met.   It  is  our  aim  to  incorporate 
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a  means  by  which  that  particular  curve  Is  singled  out  which 
also  satisfies  the  free  boundary  condition.   We  achieve  this 
by  use  of  the  variational  principle  described  in  the  next 
section. 

3. 5«  Variational  Principle  —  Free  Boundary  Problems 

We  establish  the  free  boundary  by  applying  the  method  of 
steepest  descent  to  the  following  variational  principle  associa- 
ted with  free  boundary  problems  [  5].   Among  all  admissible 

boundaries  r*.  each  of  which  defines  a  domain  D   in  which  we 

z 

can  find  a  stream  function  tl/  which  satisfies  the  Dirichlet 
conditions  on  the  boundary,  the  appropriate  free  boundary  r  is 
characterized  by  the  property  that  the  energy  M  of  the  flow  is 
a  minimum  for  a  fixed  representative  mass  (or  volume)  V.  This 
principle  can  be  formulated  as  an  isoperimetric  problem  in  the 
calculus  of  variations  in  which  we  minimize  the  functional 


Jp  =  M  -  A^ 


with  respect  to  all  admissible  free  boundary  configurations. 
As  a  particular  example,  the  appropriate  expressions  for  M  and 
V  for  the  Riabouchinsky  models  can  be  shown  to  be 


M  =  7r"^j/  [V(^-yl+"^/(l+m))]2y"^  dxdy 
D 


V  =  {2ir)'^  J^f  y"^   dxdy 


O 

where  O  represents  the  cavity  [3].   If  5n  represents  an 
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arbitrary  perturbation  of  the  free  boundary  along  its  inner 

normal  then  we  can  calculate  the  first  variation  of  J„  to  be 

F 


5jp  =  f  [iWy'^f   -A^lSn  ds  . 


Thus  we  obtain  as  a  necessary  condition  for  J_  to  be  a  minimum 


the  requirement 


(V^/y"^)2  =  a2 


which  is  the  free  boundary  condition,  and  the  Lagrange  multi- 
plier A  represents  the  constant  speed  q  along  r. 

Applying  the  method  of  steepest  descent  as  a  technique 
for  implementing  the  minimization  of  J„,  we  are  led  to  the 


formula 


m       ^2  dV 
dt 


7.^  §  =  f  [Wy")^  -^^]r,^i'^ 


where  n,  is  the  speed  of  r  in  the  direction  of  its  inner  normal, 
and  we  specify 

(6)  n^  =_y((VVy"^)^-A2) 


as  the  path  of  steepest  descent,  where  the  positive  number  y 
represents  a  convergence  factor.   This  variational  principle 
applied  to  the  other  models  leads  also  to  equation  (6)  as  the 
path  of  steepest  descent,  although  the  particular  expressions 
for  M  and  V  are  different.   Formula  (6)  is  essentially  the 
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recipe  used  for  modifying  the  free  boundary  in  [ 6  ]  and  [ 7  ] , 
with  A  chosen  as  the  average  of  IVT/z/y"^!  computed  over  the  free 
boundary.   In  the  present  context,  we  see  that  there  is  a 
mathematical  basis  for  expecting  such  a  recipe  to  converge, 
since  it  comes  from  a  variational  principle. 

J)A.    Combined  Iterative  Method 

We  use  the  two  variational  principles  we  have  described 
to  produce  the  following  basic  iterative  scheme.   In  the  fixed 
rectangle  R  we  solve 

X,  =  aAx 
(7a)  y^  =  aAy 

Tp^   =  a[^f   -m/y(VV/-Vy)] 

We  remark  that  the  time-dependent  formulation  of  the  equation 
for  Tp   can  be  motivated  by  an  application  of  the  method  of 
steepest  descent  to  Dlrlchlet's  principle  for  f   over  the 
physical  domain.   Furthermore,  the  form  taken  by  this  equation 
for  m  =  1,  the  axially  symmetric  case,  is  to  maintain  Invariance 
under  the  conformal  transformation  z(w). 

On  the  boundary  of  R  we  impose  conditions  dictated  by 
the  variational  principles.   Along  those  parts  of  ^R  which  map 
into  the  fixed  boundary  of  D,  we  have 

(7b)  x^+iy^  =  P[(%^s  +  yv^s^/(^s+y^^(^s  +  ^ys^ 
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where  the  s  and  v  subscripts  are  either  u  or  v  derivatives. 
Equation  (7h)  Is  the  tangential  shift  motivated  by  the  con- 
formal  mapping  procedure.   On  the  part  of  ^R  which  maps  Into 
the  free  boundary  we  Impose,  In  addition  to  the  tangential 
shift,  a  normal  shift  based  on  equation  (6).   We  arrange  the 
boundary  correspondences  so  that  the  base  of  the  rectangle 
corresponds  to  the  free  boundary,  hence  s  =  u  and  v  =  v,  and 
we  obtain 

(7c)   X,  +  iy ,  =  p[(x  X  +y  y  )/(x  +y  )](x  +  ly  ) 

+  7[y^"'(x^+y^)/V'^  -lA2](x^+ly^)  . 

The  second  term  on  the  right  in  (7c)  represents  an  implementa- 
tion of  equation  (6),  where  we  have  used  the  transformation 


(8)  (V^^)2  =  t^/(x^+y2)  . 


The  gradient  V  is  applied  in  the  physical  domain  along  r. 
This  transformation  is  valid  when  the  mapping  z(w)  is  conformal 
and  we  could  use  u-derlvatives  as  well  in  the  denominator 
of  (8).   The  parameter  A  is  a  suitably  chosen  constant.   The 
vector  (x  ,y  )  is  the  approximate  inner  normal  to  r. 

The  boundary  conditions  for  ^  are  the  Dirichlet  conditions 
transferred  to  the  boundary  of  R. 

Assuming  that  a  steady  state  is  reached,  we  have 


32 


Ax  =  0 
Ay  =  0 
A^  -  (m/y)VT//.  Vy  =   0 


in  R, 


(9a)  Vv+Vv=° 


2        2  2        2 

(9b)  X    +y^  =  x'^  +  y^  +H 


in  R  and  on  bR,    and 


(10)  y^"'(x^+yJ)/l^J  =  lA^  =  constant 


on  that  part  of  ^R  which  maps  into  r.   Using  (8),  we  see  that 
(10)  expresses  the  fact  that  the  free  boundary  condition  is 
satisfied.   We  have  yet  to  indicate  how  the  requirement  H  =  0 
is  achieved.   This  depends  upon  the  choice  for  A.   For  example. 
In  the  cusped  geometry  model,  (where  m  =  0)  we  are  able  to 
obtain  convergence  with  the  formula 

where  the  subscript  denotes  an  average  over  the  free  boundary. 
To  show  that  H  =  0,  we  combine  (10)  with  (9b)  to  obtain 

lA'^  =  E/f^   +  (X^  +y^)/^2  ^ 
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Upon  taking  averages  in  this  equation  and  using  the  form  for 
A  we  have  chosen,  we  get 

H  =  0  . 

For  the  hydrodynamlc  models,  the  situation  Is  consider- 
ably more  delicate  as  regards  the  choice  of  A.   We  consider 
the  necessary  refinements,  as  well  as  other  details  of  the 
method  as  applied  to  each  model.  In  Section  4. 

3.5«  Remarks  on  the  Method 

1.  We  observe  that  the  kind  of  systematic  Iterative 
procedure  we  have  described  differs  sharply  from  previous  work 
In  that  the  free  boundary  r  In  our  method  Is  modified  at  each 
time  step,  consistent  with  the  overall  formulation  of  a  time- 
dependent  problem.   In  previous  work,  a  succession  of 
Dlrlchlet  or  linear  mixed  problems  were  solved  completely  with 
a  fixed  guess  for  the  location  of  r.   The  free  boundary  was 
then  modified  by  an  appropriate  application  of  the  free  boundary 
condition  to  produce  a  better  guess. 

2.  We  have  considerable  freedom  in  our  formulation  to 
choose  the  boundary  correspondences  for  the  conformal  mapping. 
We  take  advantage  of  the  right-angled  corners  of  the  rectangle 
by  specifying  that  one  of  these  corners  map  Into  the  separation 
point  between  the  free  and  fixed  boundaries  of  the  flow  region. 
In  this  way,  both  z  and  ip   become  regular  functions  of  w  at  the 
corner  and  the  difficulty  caused  by  the  singular  behavior  of 

f   in  the  z-plane  is  removed.   The  finite  difference  analogs  of 
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the  time  dependent  equations  can  be  made  accurate  to  second 
order  in  mesh  width  h  over  the  whole  rectangle,  and  we  expect 
the  convergence  to  the  continuous  solution  to  be  0(h). 

3.  We  discuss  the  application  of  the  method  to  a 
completely  trivial  example  for  the  purpose  of  providing  a 
geometric  interpretation  of  the  ideas  we  have  introduced. 

Consider  the  problem  of  mapping  the  square  S  onto  the 
rectangle  shown  in  Figure  5»1,  such  that  the  four  corners  A, 
B,  C  and  D  of  the  square  map  into  the  corresponding  four 
corners  of  the  rectangle.   It  is  obvious  that  if  we  want  the 
mapping  to  be  conformal  then  the  value  of  a   must  be  1  and  the 
only  solution  is  the  identity 


z  =  w 


Thus  we  shall  ultimately  allow  the  boundary  BC  in  the  z-plane 
to  translate  horizontally,  governed  in  its  motion  by  the 
normal  shift  condition,  but  let  us  assume  for  now  that  this 
boundary  is  fixed  at  x  =  a.   The  affine  mapping 

y  =  V 
(11) 

X  =  au 

is  obtained  as  the  steady  state  solution  which  minimizes  the 
Dirichlet-Douglas  integral.   The  functions  x  and  y  are  harmonic 
and  we  obviously  have 
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X  X  +  y  y  =0 
u  V   ■'u  V 


2   2     2    2      2 

x^ +y^  -  x^  -y^  =  1-a  =  constant 


In  the  square  and  on  its  boundary.   We  see  that  only  if  a  =  1 
will  the  mapping  be  conformal. 

The  general  theory  of  the  Plateau  problem  asserts  that 
minimal  surfaces  in  three  dimensions  (x,y,s),  where  we  use  s 
as  the  third  coordinate  to  avoid  confusion  with  the  complex 
variable  z,    satisfy  the  parametric  equations 


in  S,    and 


Ax   =   Ay  =  As   =   0 


xx+yy+ss     =0 
u  V      ''u'^v        u  V 


2         2         2  2         2         2 

X     +y    +  s      -  X     - y    -  s     =   0 
v     "^v        V  u     -^u        u 


on  the  boundary.   In  the  present  context,  we  may  interpret  the 
solution  (11)  as  the  projection  onto  the  (x,y) -plane  of  a 
three-dimensional  minimal  surface,  whose  s-coordinate  satisfies 


(12) 

on  SS  and 


s  s  =  0 
u  V 


2   2   ,  ^2 

s  -  s   =  1-a 
u    V 


As  =  0 


in  S.   Suppose  ct  >  1.   Then  for  s  to  be  real  we  must  have, 
from  equations  (12) 
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s^=  0 


from  which  we  obtain 

s  =  (a^-l)^/^v  . 

Thus  the  three-dimensional  minimal  surface  -^b  _  ciqu^re  in  the 
coordinate  system  (x,y,s)  which  is  hinged  along  the  x-axis  and 
inclined  in  the  s-direction  an  appropriate  amount  so  that  its 
length  is  a.   Conversely,  if  a  <  1,  we  obtain 


=  (l-a2)l/2u  , 


which  gives  (as  the  minimal  surface)  an  inclined  square  of  side 
1  which  is  hinged  along  the  y-axis.  Thus  the  effect,  geometri- 
cally, of  the  normal  shift  in  our  formulas  is  to  determine  that 
minimal  surface  in  a  one-parameter  family  of  surfaces  for  which 

s  =  0  , 

in  which  case  the  mapping  will  be  conformal.   Assume  that  l/A 

2    2 

is  chosen  to  be  the  average  of  x  +y  over  BC-   If  a  >  1,  this 

normal  shift  drives  the  boundary  BC  inward,  as  we  see  by  inspec- 
tion of  equation  (3.7c).  If  a  <  1,  BC  is  moved  outward.  There- 
fore we  expect  the  combined  iterative  procedure  to  converge. 
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^.  THE  FINITE  DIFFERENCE  SCHEME 

The  basic  numerical  scheme  described  in  Section  5  is 
applied  to  the  cusped  geometry  model,  and  to  the  plane  and 
axially  symmetric  vena  contracta  and  Riabouchinsky  models 
introduced  in  Section  2.   The  cusped  geometry  model  was  chosen 
for  study  because  it  provides  a  simple  example  on  which  to  test 
the  method.   As  we  will  describe  in  Section  4.2,  the  method  is 
applied  in  this  case  as  a  systematic  iterative  procedure  for 
determining  the  appropriate  physical  domain  in  which  the  plasma 
physics  problem  is  solved.   No  attempt  is  made  to  give  special 
treatment  to  the  separation  points.   In  the  formulation  of  this 
model  we  allow  points  to  shift  back  and  forth  between  free  and 
fixed  boundaries  as  t  — »■  oo  .   As  a  result  the  over-all  method 
possesses  a  high  degree  of  stability  and  there  are  no  inherent 
problems  in  obtaining  convergence,  although  the  accuracy  is 
poor.   For  the  hydrodynamics  models  we  introduce  the  additional 
feature  of  mapping  a  corner  of  i:he  parameter  domain  rectangle 
into  the  separation  point  for  the  purpose  of  eliminating  the 
effects  of  the  singular  behavior  of  the  stream  function.   As 
we  describe  in  Section  4.3,  this  formulation  requires  procedures 
of  a  considerably  more  delicate  nature  to  produce  a  convergent 
scheme.   We  study  the  plane  models  so  that  we  can  check  our 
numerical  computations  against  the  known  analytic  solutions. 
The  axially  symmetric  models  are  calculated  to  demonstrate  that 
the  method  can  be  applied  to  examples  that  are  governed  by 
differential  equations  other  than  Laplace's  equation  in  two 
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dimensions.   We  also  obtain  for  these  cases  new  estimates 
of  certain  parameters  of  physical  importance  for  the  models 
studied.   These  are  the  contraction  coefficients  for  the  vena 
contracta  and  the  drag  coefficient  for  the  Riabouchinsky  model. 

k.l.    General  Features  of  the  Numerical  Method 

4.1.1.  Basic  Algorithm 

Finite  difference  equations  are  written  to  approximate 
the  differential  equations  for  x,  y,  and  ^  in  R  and  on  its 
boundary.   We  observe  that  since  the  basic  domain  is  a  rec- 
tangle, it  is  a  simple  task  to  produce  finite  difference  equa- 
tions which  are  accurate  to  second  order  in  mesh  size. 
Furthermore,  we  can  arrange  the  dimensions  of  the  rectangle  to 
have  a  rational  length-to-width  ratio  so  that  the  mesh  spacing 
h  is  equal  in  both  horizontal  and  vertical  directions. 

The  solution  is  obtained  by  the  following  steps: 

1.  An  initial  mapping  is  established  which  sets  up  the 
boundary  correspondences  to  conform  closely  with  the  conformal 
mapping  we  seek.   An  initial  solution  for  tp   is  also  selected. 

2.  The  interior  mesh  points  are  scanned  and  we  move  for- 
ward one  time  step  At  in  the  solution  of  equations  (3.7a)  for 
X,  y,  and  f. 

J>.    The  tangential  shift  (3.7b)  is  applied  to  those  points 
on  SR  which  map  onto  the  fixed  boundary  of  SD,  which  relocates 
their  positions  on  the  fixed  boundary  at  the  next  time  step. 
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Where  appropriate,  higher  order  normal  corrections  are  applied 
to  keep  these  points  on  the  fixed  boundary. 

h.    The  combined  tangential  and  normal  shift  given  by 
(3.7c)  is  applied  along  the  free  boundary. 

5.  A  check  for  convergence  is  performed.   We  examine  the 
interior  residuals  as  well  as  the  boundary  terms 


and 


(x  X  +  y  y  )/(x  +y  ) 

^uv  ''uv^'^u  ^u' 


y^"'(x^y^)/*2  _  lA^ 


to  see  if  the  maximum  of  their  absolute  values  lies  below  some 
specified  error  E.  If  not,  we  go  back  to  step  2.  If  so,  then 
the  iterative  process  is  terminated. 

i^-.1.2.   Interior  Equations 

The  method  of  steepest  descent  produces  the  heat  equation 
(3.7a)  to  be  solved  in  R  as  t  -►  00  .   The  formulas  which  result 
from  using  a  forward  difference  for  the  time  derivative  with 
time  step  At  =  4h  ,  where  h  is  the  mesh  width,  correspond  to 
the  familiar  Jacoby  iterative  method  for  Laplace's  equation. 
This  method  is  far  too  slow  to  be  practical.   Therefore,  we  use 
successive  point  over-relaxation,  which  yields  the  basic 
difference  equation 

where  cd  is  the  relaxation  factor  and  j  denotes  the  time  step. 
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Identical  equations  are  used  for  y  and  for  ■^  in  the  plane 
case.   The  axially  symmetric  difference  equation  for  tj/   will 
be  derived  in  Section  4.1.3.   As  t  -^  oo  we  obtain  a  finite- 
difference  solution  accurate  to  second  order  in  h. 

To  conform  to  our  point  of  view  that  the  method  specifies 
a  time-dependent  problem  which  approaches  a  steady  state,  we 
note  that  the  over-relaxation  formula  (1)  has  an  interpretation 
as  the  numerical  analog  of  the  hyperbolic  equation  [ll]. 


(2)  x,-x,+ax,  =Ax, 

^    ^  ut   vt     t      ' 


where  the  time  step  At  in  the  finite  difference  analog  is  taken 
equal  to  h.   The  constant  a,  which  we  may  interpret  as  a 
dissipation  parameter,  is  related  to  the  relaxation  factor  o) 
by  the  formula 

(5)  CD  =   2/(1  +ah)  . 

Fourier  analysis  [11]  can  be  applied  to  show  that  the  finite 
difference  solution  to  (2)  approaches  a  steady  state  according 
to  the  estimate 

(4)  |x°°   -xN   I  <  Ae-^^  , 

^  '  '  m,n   m,n'  —        ■* 

where  N  is  the  total  number  of  iterations,  A  and  A  are  numbers 
which  depend  only  on  a,  and  x    is  the  N   iterated  solution. 
This  formula  shows  that  the  number  of  iterations  required  to 
achieve  a  specified  accuracy  grows  like  l/h  provided  that  the 
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relaxation  factor  is  adjusted  in  accordance  with  formula  (3) 
to  maintain  a  constant  a.   The  convergence  rate  indicated  by 
(^)  is  in  sharp  contrast  to  that  which  is  obtained  by  the  Jacoby 
scheme.   There,  the  number  of  Iterations  grows  like  l/h   [1  ]. 
In  the  present  context  the  improvement  provided  by  over- 
relaxation  is  due  to  the  fact,  indicated  above,  that  we  can 
take  At  =  h  in  the  numerical  analog  of  equation  (2). 

4,1.5.  The  Axially  Symmetric  f   Equation 

In  the  axially  symmetric  models  the  stream  function 
satisfies  the  equation 

(5)  A^^  -  ^y/y  =  0 

in  the  z -plane.   Invariance  under  a  conformal  transformation 
z  =  z(w)  is  obtained  when  this  equation  is  written  in  the  form 

(6)  At//  -  (Vy  .  VV/)/y  =  0  . 

Thus  we  may  regard  ip   and  y  in  equation  (6)  as  functions  of  w, 
with  the  operators  A  and  V  acting  in  the  w-plane.   We  seek  ^  as 
a  function  of  w  over  the  rectangle  R  in  our  parametric  formula- 
tion.  As  shown  below,  we  can  write  finite  difference  equations 
for  (6),  which  are  accurate  to  second  order  in  mesh  size  and 
which  are  suitable  for  solution  by  over-relaxation.   However, 
for  both  hydrodynamic  models,  the  conformal  mapping  we  seek 
between  R  and  the  physical  plane  has  a  pole  at  a  corner  of  R. 
In  the  neighborhood  of  this  pole  the  truncation  error  in  the 
finite  difference  approximation  to  (6)  becomes  unbounded. 
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We  can  remove  this  difficulty  by  use  of  the  Kelvin 
inversion  [12] 


where 


and 


^*(X,Y)  =  iz|^(x,y)  , 
Z  =  X  +  iY  =  l/z 


|Z|^  =  X^+Y^  =  l/(x^+y^)  . 


If  Tjj   satisfies  equation  (5),  then  ip*   satisfies 


A^t*  -  ^y/^  =   0    ,  .. 


or,  in  the  w-plane 

( 7 )  A^*  -  VY  •  V^*/Y  =  0  . 

The  transformation 

Z  =  l/z 

maps  the  neighborhood  of  infinity  in  the  z-plane  onto  a  neigh- 
borhood of  the  origin  in  the  Z-plane,  and  equation  (?)  can 
therefore  be  used  instead  of  equation  (6)  in  the  vicinity  of 
the  pole. 

Finite  difference  equations  for  either  (6)  or  (7)  which 
are  accurate  to  second  order  in  h  can  be  obtained  following 
the  procedure  in  [7  ].   We  write  the  differential  equation  in 
divergence  form 
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yt|i  (Vy>  +l7(V=''i  =  °  • 


Let  the  points  in  a  grid  of  uniform  mesh  of  width  h  be  numbered 
as  shown  in  Figure  4.1.   A  central  difference  on  the  outer 
u-derivative  gives  the  expression 

^o^^V^^a  -  (Vy^b^/^^^/^^ 

for  the  derivative  at  the  point  0.   If  we  use  central  differ- 
ences on  Tp     at  the  points  a  and  b  and  simple  averages  for  y 
at  a  and  b  relative  to  the  points  0,  1  and  3,  0  respectively, 
we  obtain  the  formula 


(8)   A  =  ay^h^ 


{^i-^o)/(yi+yo)  -  i^^--p^)/{Yj,+y^) 


as  an  approximation  to  y  -^  {'^  ./v)    at  the  point  0.   To  see  that 
it  is  accurate  to  second  order  in  h,  we  note  that  A  considered 
as  a  function  of  h  has  the  expansion 

2 

A  =  Aq  +A^h  + Agh  +  .  ..  . 

Since  the  quantities  f   and  y  satisfy 


y^(h)  =  y-L(-h)  , 


we  conclude  by  inspection  of  the  terms  which  appear  in  equation 
(8)  that  A  is  an  even  function  of  h.   In  particular. 
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A^  =  0  . 


Direct  computation  using  Taylor  series  expansions  shows  that 

^o  =  ^^uu  -  ^uV^^o 

p 

and  so  A  is  a  finite  difference  approximation  accurate  to  0(h  ) 

Taking  also  the  v-derivative  into  account,  the  finite 
difference  equation  for  ijj   at  the  point  0  is 

(2y  /h^)^  {f^-f^)/{Y^+y^)    =  0  . 
1=1 


We  solve  this  equation  for  ij/   ,    to  obtain  the  form 


1=1 
where 


4 


1    '"  1  ''o 
and 


a,  =  (y,+y^)  "^  ,        i  =  1,2,3,^, 


4 


which  exhibits  the  strict  diagonal  dominance  of  the  matrix 
for  f.      This  matrix  is  also  symmetric  and  therefore  the  system 
of  equations  is  suitable  for  solution  by  successive  over- 
relaxation.   There  are  no  difficulties  associated  with  the 
circumstance  of  solving  near  y  =  0. 
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4.1.4.  Boundary  Conditions 

All  space  derivatives  at  the  boundary  of  R  are  calculated 
to  second  order  In  h,  using  a  central  difference  for  the 
tangential  derivatives  and  a  3-polnt  formula  for  the  normal 
derivatives.   These  formulas  are  given  below  relative  to 
Figure  4.2,  which  shows  appropriately  labelled  points  In  a 
grid  of  mesh  width  h. 


^9a)        U 


(9^)       i 


=  -{kf^-3f^  +  f^)/2h  +   O(h^)  , 
=  {t^  -f5)/2h  +  O(h^)  . 


The  quantity  f  stands  for  x,  y  or  f.      These  derivatives  are  used 
to  calculate  the  right-hand  sides  of  the  boundary  equations 
(3.7b)  and  (3' 7c).   For  the  time  derivatives  In  these  boundary 
equations  we  use  simple  forward  differences  and  we  take  the 
time  step  At  equal  to  h.   Thus  we  obtain  difference  equations 
at  the  boundary  of  the  form 

where  z  =  x  + ly,  the  superscript  J  denotes  the  time  step,  z^ 
and  z^  are  the  numerically  computed  tangential  and  normal 
derivatives  at  time  step  J,  T*^  and  N^  are  the  numerically 
calculated  tangential  and  normal  terms  In  square  brackets  in 
equations  (3.7b)  and  (3»7c)  which  approach  zero  as  j  -»■  cd  ,  and 
P  and  7  are  the  tangential  and  normal  convergence  factors. 
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These  factors,  taken  sufficiently  small  to  produce  a  convergent 
scheme  for  a  particular  mesh  size,  are  held  fixed  under  mesh 
refinements  so  that  the  finite  difference  equations  at  the 
boundary  remain  the  numerical  analogs  of  the  corresponding 
differential  equations. 

4,1. 5^  Convergence  Rate  of  the  Iterative  Scheme 

The  convergence  factors  P  and  y   in  the  boundary  equations 
and  the  dissipation  parameter  a  in  the  interior  equations  must 
be  chosen  to  produce  a  convergent  scheme.   The  theory  of  the 
method  of  steepest  descent  indicates  that  these  quantities  have 
to  be  taken  sufficiently  small  in  order  to  avoid  over-shooting 
the  minimum  point  on  the  multi -dimensional  surface  described 
by  the  discrete  formulations  of  the  Dirichlet-Douglas  func- 
tional 3«1  and  the  free  boundary  functional  J_  of  Section  3.3. 
Assuming  that  these  parameters  are  appropriately  chosen  to 
produce  convergence,  linearized  theory  suggests  that  for 
sufficiently  small  mesh  size  the  estimates 


|x°°   -xN   I  <  Ae-^^ 
m,n   m,n'  — 


|y^  -y^  I  <  Ae-^^ 
m,  n   m,  n'  — 


^°°  -r      I  <  Ae 
^m,n  ^m,n'  — 


are  valid  for  the  rate  of  convergence  to  the  steady  state 
finite  difference  solution,  just  as  estimate  (4)  holds  for 
equation  (2).   Here  A  and  A  are  constants  independent  of  h  and 
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A  is  related  to  the  lowest  eigenvalue  of  the  linearized 
problem.   If  a,  p,  and  7  are  held  fixed  under  mesh  refinement, 
then  we  expect  the  number  of  iterations  N  to  grow  like  l/h. 
We  should  observe,  in  this  connection,  that  this  rate  is  the 
fastest  that  can  be  expected  when  the  boundary  is  shifted  using 
a  forward  difference  for  the  time  step.   In  particular,  there 
is  no  advantage  to  using  more  sophisticated  schemes  (such  as 
alternating-direction  methods)  on  the  interior  equations, 
because  the  boundary  shifts  will  still  keep  the  overall  con- 
vergence rate  at  0(l/h).   Furthermore,  the  computation  time  T 
increases  like  l/h-^,  because  the  number  of  mesh  points  grows 
like  l/h  ,  and  the  total  computing  time  is  givfen  approximately 
by 

T  =  kMN  , 

where  k  is  the  time  to  compute  at  a  single  mesh  point,  M  is 
the  total  number  of  mesh  points  and  N  is  the  total  number  of 
iterations  required  for  convergence. 

4.2.  Cusped  Geometry  Model 

The  numerical  method  was  applied  first  to  the  cusped 
geometry  problem  described  in  Section  2.1  because  it  provides 
a  relatively  simple  example  on  which  to  check  the  validity  of 
the  scheme.   In  this  model  no  attempt  is  made  to  map  corners 
of  the  rectangle  into  the  separation  points.   Also,  the  f 
equation  does  not  enter  into  the  calculations.   Thus  we  are 
dealing  essentially  with  the  problem  of  obtaining  a  suitable 
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periodic  conformal  mapping  of  the  rectangle  onto  an  annulus 
whose  outer  boundary  Is  constrained  In  a  special  way. 

4.2.1.  Formulation  of  the  Equations 

Considerations  of  symmetry  permit  the  problem  to  be 
formulated  as  a  mapping  from  a  rectangle  R  onto  the  physical 
domain  D^  In  the  first  quadrant  of  the  z -plane  as  shown  In 
Figure  4.3.   The  base  AB  of  the  rectangle  Is  required  to  map 
onto  the  outer  boundary  T  + Tp  of  the  annulus,  where  r  Is  the 
free  boundary  and  Tp  Is  the  fixed  boundary  consisting  of 
portions  of  the  walls  x  =  i  and  y  =  i.   The  parameter  H   is 
arbitrary,  and  defines  the  basic  periodic  structure  of  the  model 
The  points  S^  and  Sg  denote  the  separation  points  between  r  and 
Pp.   Their  location  on  the  walls  is  unknown  and  must  be  found 
as  part  of  the  solution.   Our  problem  is  to  determine  the  con- 
formal  mapping  of  R  onto  D  which  has  the  boundary  correspondence 
indicated  in  Figure  4.3,  and  for  which  r  is  located  so  that  the 
free  boundary  condition  holds  along  it.   We  do  this  by  solving 
the  numerical  analog  of  the  following  system  of  time-dependent 
equations  which  result  from  the  general  theory: 

(lOa)  X,  =  oAx  ,   y  =  aAy 

In  R.   In  practice  these  equations  are  modified  by  the  use  of 
successive  over-relaxation,  as  described  in  Section  4.1.2.   On 
the  boundary  of  R  we  obtain  from  the  general  theory 
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(10b)     x,  +  iy,  =  -p[(Vv+yuyv)/^^u+yu5nx^+iyj 


along  CD,  where  the  convergence  factor  p  is  positive,  and 


(10c)  x^+iy^  =  Pt^Vv  +  Vv^/^^u+yu^lf^u+lyu^ 


+  7[x2+y2   -    {x^  +  y2)J(x^  +  iy^) 


along  AB.   On  the  segments  AS-,  and  S-,B  which  are  supposed  to 
map  onto  the  walls  r„  we  effectively  set  7  equal  to  zero  so 
that  these  points  are  subjected  to  a  tangential  shift  only. 
The  segment  S-,  Sp  maps  onto  the  free  boundary  r  and  here  both 
normal  and  tangential  shifts  are  applied  with  positive  conver- 
gence factors  P  and  7.   We  observe  that  the  location  of  the 
images  S,  and  Sp  must  be  found,  that  is,  it  is  not  known  which 
points  on  AB  map  onto  r  and  which  map  onto  r^.   The  correct 
apportionment  is  determined  as  the  problem  converges  to  a 
steady  state.   Provision  must  be  made  which  allows  points  to 
move  to  either  r  or  r  as  t  -*  00  .   We  describe  the  technique  for 
doing  this  in  Section  h.2.J>.      The  remaining  boundary  conditions 
are 


(lOd)  ^u  =  °  '    y  ^  ° 


on  AD,  and 


(lOe)  X  =  0  ,   y^  =  0 


on  EC.   These  conditions  are  Imposed  by  symmetry.   They  are 
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obvious  consequences  of  the  Schwarz  reflection  principle,  but 
can  also  be  derived  from  our  tangential  boundary  conditions. 
For  example,  on  AD,  since  we  have 

y  =  O  , 
we  also  have 

so  that  the  tangential  equation  (5.7b)  reduces  in  this  case  to 


X.  =  Px 
t   "^  u 


y^  =  0   (y  =  0)  . 


Since  we  expect  x,  -♦  0  In  the  steady  state,  we  simply  set  x 
equal  to  zero  at  each  iteration  by  an  appropriate  reflection 
in  the  finite  difference  equations  for  (lOa). 

The  stream  function  does  not  appear  in  our  formulation. 
As  shown  in  Section  2.1,  its  solution  to  the  Dirichlet  problem 
in  the  rectangle  is  simply 

^{u,v)  ="v/b  . 

The  stream  function  enters  into  the  conformal  mapping  problem 
only  through  the  normal  shift  term  along  the  free  boundary. 
Since  in  this  model  we  have 

fy      =   1/b    , 

the  constant  b  can  be  absorbed  into  the  convergence  factor  7  so 
that  equation  (10c)  contains  no  trace  of  f. 
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4.2.2.  Initial  Mapping 

The  essential  features  required  for  the  Initial  mapping 
are  that  It  be  univalent  and  that  It  possess  the  correct 
boundary  correspondences  to  the  extent  possible.   The  require- 
ment for  unlvalence  Is  necessary  because  there  Is  no  explicit 
mechanism  In  the  numerical  scheme  to  prevent  overlapping  of 
points,  particularly  those  on  the  boundaries.   The  requirement 
for  the  appropriate  Initial  boundary  correspondences  means  that 
the  Initial  mapping  takes  the  line  CD  of  the  rectangle  onto  the 
arc  CD  of  the  physical  plane,  and  takes  the  vertical  sides  onto 
the  corresponding  parts  of  the  physical  plane  as  Indicated  In 
Figure  4.3.   The  Image  of  the  base  AB  of  the  rectangle  cannot 
be  specified  as  completely.   It  is  precisely  Its  location  which 
is  to  be  found  as  t  -♦  oo  . 

In  practice  we  chose  the  rectangle  dimension  a  =  7r/2 
and  b  equal  to  a  rational  multiple  of  a  so  that  mesh  spacing 
in  horizontal  and  vertical  directions  can  be  taken  equal.   The 
inner  radius  r  was  chosen  to  be 


-b 
r  =  e 


With  this  scaling  the  analytic  function 

(11)  z  =  e^^ 

maps  R  onto  the  part  of  the  annulus  r  <    \z\   ^1   which  lies  in 
the  first  quadrant,  so  that  for  i  >  1,  (11)  is  the  trivial 
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solution  to  the  problem.   For  non-trivial  solutions  we  must 
choose  H   in  the  range 

V  <   H   <   1    , 

This  requires  the  physical  domain  D   to  distort  by  pushing  out- 
ward,  establishing  the  boundaries  r^  along  the  walls  x  =  ^, 
y  =  i,  and  yielding  a  conformal  mapping  with  the  same  conformal 
type  as  the  annulus. 

A  possible  initial  mapping  which  satisfies  the  require- 
ments given  above  is 

z  =  [(^-r2)e"''+r^(l-i)e'']/(l-r^)e^'^ 

which  takes  R  onto  the  part  of  the  annulus  v  <_    |z|  <  I   which 
lies  in  the  first  quadrant.   Its  real  and  imaginary  parts  are 
harmonic  and  the  steady  state  tangential  condition 


X  X  +  y  y  =0 
u  v  *^u''v 


is  satisfied.   For  the  normal  shift,  we  have,  initially, 

(12)    4+yv-f^u+yu^o  =  -2r^(l-^)(^-r^)/(l-r'^)^  <   0  , 

and  since  this  quantity  is  negative,  the  outer  circle  |z|  =  I, 
which  is  the  initial  guess  of  r,  will  be  pushed  outward  toward 
the  walls  x  =  i,  y  =  i.   The  technique  for  dealing  with  the 
constraint  that  points  must  not  move  beyond  these  walls  is 
described  in  the  next  section. 
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^.2.5.  Boundary  Conditions 

When  steady  state  is  achieved  in  equations  (10),  the 
free  boundary  condition 


2,2   ,2,2, 
X  +y  -  fx  +y  J   =  0 


is  satisfied  on  r,  and  the  mapping  is  conformal  in  consequence 
of  the  general  theory  described  in  Section  3»^»   For  the  way 
in  which  we  have  formulated  the  cusped  geometry  model,  pro- 
vision must  be  made  that  the  proper  apportionment  of  points 
between  the  fixed  boundary  Tt-,  and  the  free  boundary  r  evolves 
as  t  — ^  00  .   To  describe  our  method  for  achieving  a  stable  con- 
figuration of  the  outer  boundary  T  +  Tt:,,  ^^  define  S  to  be  the 
domain  bounded  by  the  coordinate  axes,  the  lines  x  =  ^,  y  =  i 
and  the  circle  |z|  =  r.   As  we  have  observed  in  the  previous 
section,  if  i  >■  1  then  the  solution  we  seek  reduces  to  the 
trivial  conformal  mapping 

iw 
z  =  e 

for  the  scaling  we  have  chosen.   If  ^  <  1  then  it  is  clear 
that  the  physical  domain  D  must  fill  out  an  appropriate  por- 
tion  of  S  in  order  to  be  the  image  by  a  conformal  mapping  of 
R,  which  has  the  same  conformal  type  as  the  above  trivial  solu- 
tion.  For  this  conformal  mapping  a  certain  number  of  points 
along  AB  will  have  their  images  on  the  walls  x  =  £,    y  =  £, 
forming  the  fixed  boundary  r„,  and  the  rest  will  map  onto  r, 
which  will  be  located  in  the  Interior  of  S.   We  describe  two 
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possible  extreme  configurations  which  demonstrate  that  the 
normal  shift  controls  the  convergence  to  the  mapping  with  the 
correct  conformal  type.   Consider  the  Initial  mapping  given 
in  the  last  section.   Here  Tp  consists  only  of  the  endpoints 
^  and  \.l,    where  I   <   \.      As  we  see  by  inspection  of  equation 
(12),  the  normal  shift  term  tends  to  push  the  initial  free 
boundary  outward.   This  motion  results  in  points  near  the  ends 
of  the  outer  boundary  becoming  attached  to  walls  so  that  they 
become  part  of  r„.   On  the  other  hand,  consider  a  mapping  which 
fills  out  a  substantial  portion  of  S  so  that  r  consists  of  a 
sharply  rounded  segment  near  the  corner  i+i^  of  S,  and  suppose 
that  all  steady  state  equations  in  R  and  on  its  boundary  are 
satisfied  with  the  exception  of  the  normal  shift  term.   We  use 
a  local  argument  which  suggests  that  the  normal  shift  will 
push  points  inward,  away  from  the  corner.   The  derivatives  with 

respect  to  u  and  v  are  high  along  the  rounded  segment; 

2   2 

specifically,  they  are  larger  than  the  average  (x^+yu)o  ^^^^"^ 

over  r,  so  that  the  normal  shift  drives  points  inward,  and  in 
particular  it  will  move  points  away  from  the  walls  so  that  they 
become  part  of  r. 

Between  these  two  extreme  mappings  it  is  plausible  that 
a  stable  configuration  exists  with  the  proper  conformal  type. 
We  build  into  the  formula  (10c)  by  which  the  boundary  points 
along  AB  are  shifted  a  mechanism  which  senses  the  current  state 
of  the  outer  boundary  and  adjusts  points  according  to  rules 
which  are  suggested  by  the  behavior  of  the  extreme  cases 
examined  above.   At  time  t,  points  along  AB  have  their  images 
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distributed  along  the  current  boundaries  r  and  r^^.   The 

2    2 
average  (x  +y  )   is  computed  over  those  points  which  map 

onto  r.   Using  this  value,  normal  and  tangential  shifts  are 
computed  for  the  image  of  every  point  on  AB.   If  the  shifted 
Image  lies  within  S,  then  it  is  part  of  the  free  boundary  r. 
In  particular,  a  point  previously  on  T-p  can  become  part  of  r 
in  this  way,  if  the  effect  of  the  normal  shift  is  to  drive  it 
away  from  the  wall  into  S.   If  the  shifted  point  lies  outside 
of  S,  then  the  normal  shift  term  is  dropped,  and  the  shift  is 
modified  so  that  the  image  is  pasted  to  the  wall  at  a  location 
determined  by  the  tangential  shift  alone.   Such  a  point  is 
part  of  the  fixed  boundary  ?„.   A  point  previously  on  r  can 
become  part  of  r^  In  this  way,  if  the  normal  shift  term  tends 
to  drive  it  outside  of  S.   We  observe  that  this  recipe  pro- 
vides a  technique  for  attaching  points  to,  or  detaching  points 
from,  the  walls  x  =  ^,  y  =  ^,  in  order  to  achieve  the  proper 
conformal  type.   The  locations  of  the  separation  points  S^  and 
So  between  r  and  T-n,  are  determined  as  the  solution  approaches 

d.  I* 

a  steady  state.   We  found  this  scheme  to  be  highly  stable. 
However,  there  are  large  truncation  errors  in  the  vicinity  of 
the  separation  points  S-,  and  Sp  because  no  provision  has  been 
made  to  deal  with  the  singular  behavior  there.   Moreover,  the 
discrete  images  which  define  the  endpoints  of  r  are  really  only 
approximations,  presumably  lower  bounds,  to  the  locations  of 
the  actual  separation  points  of  the  continuous  model. 
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The  algorithm  for  treatment  of  these  boundary  points 

at  time  t  Is  given  below.   At  this-  stage  of  the  computation 

we  have  available  all  derivatives  necessary  to  apply  equation 

(lOc).   We  also  have  tags  on  each  boundary  point  to  determine 

which  are  on  r  and  which  are  on  Tr^* 

2    2 
1]   Compute  the  average  (x  +y  )   for  those  points  currently 

tagged  as  members  of  r. 

2]   For  each  point  z  on  AB: 

a)  Compute  the  point  z'  given  by  the  finite  difference 
analog  of  equation  (lOc), 


z'  =  Z+Zj^  +  Z^, 


where  z„  and  z^  are  the  total  normal  and  tangential  shifts. 
NT  ^ 

b)  Compare  Re  (z')  and  Im  (z')  with  i: 

1)   If  Re  (z»)  ^  -«  and  Im  (z')  ^  -« 
then  tag  z '  as  a  point  on  r. 

11)   If  Re  (z>  )  1  ^ 
then  modify  z '  by  the  formula 


z'  =  ^+1  Im  {z+z„) 


and  tag  z'  as  a  point  on  r^. 

if 

111)   If  Im  (z')  ^  -2 
then  modify  z '  by  the  formula 
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z'  =  Re  (z  +  z„) +  ±£ 

and  tag  z'  as  a  point  on  r„. 

We  remark  that  the  above  considerations  are  in  sharp 
contrast  to  the  situation  which  prevails  if  the  value  of  the 
constraint  i  is  taken  greater  than  1  so  that  the  solution  is 
the  trivial  one.   In  this  case  the  entire  outer  boundary  of 
the  annulus  is  the  free  boundary.   In  terms  of  the  above 
algorithm,  the  cases  2b  ii)  and  2b  iii)  never  occur.   The 
separation  points  coincide  with  the  corners  A  and  B,  and  the 
model  is  quite  similar  to  the  rectangle  mapping  discussed  in 
Section  3'5« 

The  boundary  CD  of  the  rectangle  maps  onto  the  circle  O 
of  radius  r  in  the  physical  plane.   The  tangential  shift  (10b) 
is  applied  to  each  point  z  on  this  boundary  and  yields  a  new 
point  z-|  which  is  not  necessarily  on  Q  because  of  the  finite 
difference  approximations.   A  second  order  normal  shift  of  the 
form 

z'  ='r  z-^/lz^I 

is  applied  to  obtain  the  new  point  z'  on  O. 

4.3.  Hydrodynamics  Models 

The  application  of  the  numerical  method  to  the  solution 
of  the  hydrodynamics  models  includes  an  important  additional 
feature,  not  present  in  the  cusped  geometry  model.   For  these 
examples,  the  conformal  mapping  is  set  up  so  that  the  separation 
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point  z  in  the  physical  plane  Is  the  Image  of  one  of  the 

corners  w  of  the  rectangle.   This  means  that  the  conformal 
o  ° 

mapping  has  the  expansion 

p 

z  =  z+a(w-w)  +.,, 


about  the  corner.   Since  the  theory  shows  that  the  stream 

1/2 
function  is  regular  in  powers  of  (z-z  )  '      near  the  separation 

point,  we  obtain  the  expansion  for  ^(w)  at  the  corner  in  the 

form 


^  =  a^  +Im  [a^Cw-w^)  ^z.^{^^-\i^)^  +...] 


in  the  plane  case.   (The  coefficient  a^  is  actually  zero  since 

the  speed  at  z  is  finite.)  A  similar  result  holds  in  the 

axlally  symmetric  case.   The  effect  of  this  feature  of  the 

conformal  mapping  is  to  make  z  and  ^  regular  functions  of  w 

about  w  ,  so  that  the  over-all  finite  difference  method  is 
o 

accurate  to  second  order  in  mesh  width  h.   Thus  we  expect  to 
be  able  to  approximate  the  continuous  solution  as  closely  as 
we  like  by  taking  a  sufficiently  fine  mesh.   On  the  other  hand, 
the  introduction  of  this  procedure  makes  the  problem  of  pro- 
ducing a  numerically  convergent  scheme  considerably  more  deli- 
cate by  comparison  to  the  cusped  geometry  model  which,  though 
crude,  is  highly  stable.   We  indicate  the  difficulties  and  our 
methods  for  dealing  with  them  in  the  subsequent  sections. 

We  also  solve  axlally  symmetric  hydrodynamics  models  to 
show  that  our  method  can  be  applied  more  generally.   We  find 
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that  there  are  no  significant  additional  considerations  involved 
in  making  the  scheme  convergent  for  this  case,  beyond  obvious 
programming  changes.   In  the  following  sections  we  give  the 
details  of  the  formulation  for  the  vena  contracta  model.   The 
application  to  the  Riabouchinsky  model  is  essentially  the  same. 

4.3.1.  Formulation  of  the  Equations 

We  set  up  the  problem  for  the  truncated  model  described 
in  Section  2.3.   We  seek  the  conformal  mapping  from  the  rec- 
tangle R  onto  the  physical  domain  D  with  the  boundary  corre- 
spondences  shown  in  Figure  4.4. 

In  accordance  with  our  numerical  procedure,  we  solve  the 
numerical  analogs  of  the  following  equations: 

X,  =.  aAx 
{13a)         y^  =  o'Ay 

V't  =  a(A^  -  {m/y)Vy.V^ 

in  R.   In  practice,  these  equations  are  modified  by  the  applica- 
tion of  over-relaxation,  as  described  in  Section  4.1.2.   On  the 
boundary  CM  we  have 

X  =  ^ 


(13b)  ^u  =  ° 


^  =  0 
^u 


and  on  BC  we  have 
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(13c) 


^v 

— 

0 

y 

= 

0 

f 

— 

0 

These  conditions  result  from  considerations  of  symmetry.   Be- 
cause of  the  way  we  treat  the  point  at  infinity,  the  first  of 
equations  (l^b)  is  actually  different  from  the  one  shown.   We 
will  describe  the  change  in  detail  in  Section  4.2.4.   On  MA +  AB 
the  boundary  condition  for  the  stream  function  is 

^  =  1  , 

in  accordance  with  the  normalization  introduced  in  Section  2.J>. 
For  the  conformal  map,  the  tangential  condition  reduces  to 


(13d) 


X  =  0 


yt  =  -p^u 


along  AB,  using  the  fact  that  x^  is  identically  zero.   On  the 
boundary  MA,  which  is  supposed  to  map  onto  the  free  boundary  r 
in  the  physical  plane,  we  have 


(13e) 


\  +  i^t 


=  p[(x^x^  +  y^y^)/(x2  +  y2)](x^+ly^) 


+  7 


-y^"^(x^y^)/< 


lA 


2" 


(^v  +  ^^v^ 


The  appropriate  choice  for  \   is  described  in  Section  4.3.5.   In 
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the  formulation  given  by  equations  (13)*  we  observe  that  the 
ends  of  the  free  boundary  A  and  M  in  the  physical  plane  are 
permitted  to  slip  in  the  vertical  direction.   The  point  M  is 
actually  shifted  in  response  to  the  normal  term  in  (15e)  only, 
since  the  tangential  shift  is  always  zero  at  M,  by  symmetry. 
The  separation  point  A  receives  a  tangential  shift  according 
to  a  special  formula.   The  adjustment  of  the  separation  point 
is  quite  critical  because  its  location  determines  the  conformal 
type  of  the  mapping,  which  is  a  global  characteristic.   We 
discuss  this  aspect  of  the  problem  in  Section  4.3.5« 

4.3.2.  Initial  Mapping 

As  with  the  cusped  geometry  model  the  important  require- 
ments for  the  initial  mapping  are  that  it  be  univalent  so  that 
there  is  no  initial  overlapping,  and  that  it  produce  the  proper 
bo\indary  correspondence  to  the  extent  possible.   In  practice, 
the  initial  data  was  set  up  in  either  of  two  ways. 

1)  We  utilized  the  solution  to  the  plane  infinite  jet 
model  given  by  equation  (2.9)  coupled  with  the  analytic  function 

(14)    C(w)  =  I/tt  log  [-(1+cos  h(w-a)TrA)/2]  , 

which  maps  the  rectangle 

R:     0  <   Re  (w)  <   a 

0  <  Im  (w)  <  b  , 

onto  the  semi-infinite  strip 
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D^:     -00  <  Re  (O  ±   C(0,v) 
0  <  Im  (C)  <  1  . 


Examination  of  the  boundary  correspondence  shows  that  the 
composite  function  given  by  (2.9)  and  (l4), 

z  =  z(C(w))  , 

maps  R  onto  a  portion  of  the  flow  region  of  the  infinite  Jet 
which  excludes  only  the  tail  of  the  jet  as  x  -►  co  .   The  effect 
of  omitting  this  piece  is  to  truncate  the  infinite  jet  down- 
stream from  the  orifice  at  a  location  x  =  i  which  goes  to 
Infinity  linearly  with  a.   The  composite  function  and  equation 
(l4)  provide  initial  data  for  f   and  z  in  terms  of  w  which  satisfy 
all  steady  state  equations  for  the  plane  model  except  where 
the  jet  is  truncated.   We  set  the  initial  value  of  the  jet 
length  £   by  putting 

x(0,v)  =  Re  (z(0,v))  =   £    . 

2)  The  second  method  of  setting  up  initial  data  was  used 
when  the  mesh  was  refined  by  a  factor  of  two.   Given  the  numeri- 
cal solution  for  a  mesh  h,  we  interpolated  values  at  the  new 
mesh  points  to  obtain  initial  data.   This  method  was  especially 
useful  for  the  axially  symmetric  model  where  the  truncated 
plane  jet  was  actually  poor  Initial  data  for  refined  meshes. 
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k.J>.J>.    Conditions  at  Infinity 

The  conformal  mapping  we  seek  takes  the  corner  B  of  the 
rectangle  Into  Infinity  of  the  z-plane,  as  Indicated  In 
Figure  4.4.   Analysis  of  the  boundary  correspondences  shows 
that  z(w)  has  a  simple  pole  at  the  point  w„  =  a+lb.   Also, 
the  stream  function  ^  Is  discontinuous  at  w„,  with  the  value 

ID 

zero  along  BC  and  unity  along  AB.   We  deal  with  this  behavior 
by  Introducing  an  Interface  In  R,  shown  by  the  dotted  lines 
In  Figure  4.4,  which  defines  a  triangular  subdomaln  R, .   In  R, 
the  equations  are  modified  by  subtracting  off  the  singular 
behavior  of  z  and  f.      Specifically  we  define  the  functions 


(15a)  ^R^^^  "  ^^^^  "  P/^w-Wg) 


where  p  is  the  pole  strength,  and 


(15b)  ^r(w)  =  V'(w)  -f^ 


where 


^^  =  2(1-  e/TT) 


and 

^^  =  1  +  cos  Q 

As  usual,  m  =  0  corresponds  to  plane  flow  and  m  =  1  corresponds 

to  axlally  symmetric  flow,  and  9  is  the  angle  in  the  polar 

representation  of  a  point  in  D  .   The  functions  f     describe 

z  m 

the  asymptotic  behavior  of  f   as  |z|  -+ 00  in  the  second  quadrant. 

In  R,  we  solve  numerically  for  ^_  and  z„.   The  equations 

for  z„  and  il/-.   are  the  same  as  for  z  and  tj/,   because  the  functions 
n       n 
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subtracted  off  are  themselves  solutions  to  the  corresponding 

interior  differential  equations.   The  boundary  conditions 

for  z„  and  f^   take  the  form 
n       ri. 

^  =  0 


on  the  vertical  boundary  of  R-,  ,  and 


<^r'u  =  ° 


yj,  =  o 


*R  =  ° 


on  the  horizontal  boundary  of  R^.   The  boundary  conditions  on 
the  interface  are  determined  by  an  overlapping  of  the  boundary 
between  R  and  R. .   This  Interface  consists  of  two  fixed  parallel 
lines  i-,  and  H^   which  run  diagonally  through  the  mesh  as  shown 
in  Figure  4.5.   An  iteration  of  equations  (l^a)  at  a  given  time 
step  is  carried  out  in  the  following  sequence: 

1)  The  values  of  z  and  ip   on  £^   are  used  as  boundary 
conditions  in  an  Iteration  of  the  mesh  points  of  R-R-^^. 

2)  Transformations  (15)  are  applied  to  z  and  f   along  £-^ 
and  .2p  to  compute  z„  and  ^„. 
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5)  The  values  of  z„  and  V'r;  along  i,  serve  as  the  boundary 
conditions  in  an  iterative  sweep  through  the  mesh  points  of  R-,  . 

4)  Transformations  (15)  are  applied  along  £^   and  £^   to 
replace  z-^   and  f^   by  z  and  f   again. 

The  process  is  repeated  at  each  iterative  cycle.   The 
values  of  z  and  f   are  altered  along  ^,  by  step  1  and  along  Z^ 
by  step  3.   As  t-^  oo  we  approach  a  steady  state. 

In  the  axially  symmetric  case  the  additional  transforma- 
tion 

described  in  Section  4.1.3  is  used  in  R  in  order  to  keep  the 
truncation  error  0(h  ).   Thus  we  actually  solve  numerically 
for  ^*  in  R^. 

The  pole  strength  p  is  an  arbitrary  parameter  of  the 
problem  when  infinity  is  treated  in  this  way.   In  practice  we 
chose  its  value  to  be  8b/Tr{Tr+2),  which  is  the  pole  strength  in 
the  analytic  solution  for  the  plane  infinite  jet.   This  value 
is  used  for  both  plane  and  axially  symmetric  models. 

4.3.4.  Shift  of  the  Line  of  Symmetry  £ 

The  strength  of  the  pole  becomes  an  input  parameter  as  a 
result  of  our  method  of  dealing  with  infinity.   Therefore  the 
problem  as  formulated  by  equations  (13)  is  over-determined. 
We  compensate  by  permitting  the  entire  line  of  symmetry  x  =  ^ 
to  shift  horizontally,  which  allows  the  jet  length  to  be 
variable.   This  shift  is  imposed  according  to  the  formula 
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(13f)  ^t  =-   ^^\  +  ^l   -   ^u-^u^o  ' 

where  5  is  a  positive  convergence  factor  and  the  subscript  on 
the  right  indicates  an  appropriate  average  over  R.   Equation 
(15f),  which  we  call  the  ^-shift,  replaces  the  first  of  equa- 
tions (15b)  as  the  boundary  condition  on  the  vertical  side  CM 
of  the  rectangle.   Our  reason  for  expecting  a  recipe  such  as 
(13f)  to  converge  is  similar  to  the  analysis  given  in  Section 
5.5  for  rectangular  domains.   Suppose  it  is  possible  to  achieve 
a  steady  state  solution  of  the  interior  equations  and  the 
tangential  boundary  conditions.   The  theory  asserts  that  the 
quantity  x^+y^  -  x^  -  y^,  which  we  denote  by  H,  is  identically 
constant.   Suppose  that  the  value  of  &   in  this  solution  is  too 
large.   Then  we  may  assume  that  there  exists  a  rectangle  r', 
whose  length-to-width  ratio  is  larger  than  that  of  R,  which 
our  solution  maps  conformally  onto  the  flow  region.   The  mapping 
from  R  to  r'  is  an  affine  transformation  which  stretches  points 
in  the  u-direction.   Therefore  we  observe  geometrically  that 
points  in  the  tail  of  the  Jet  are  stretched  in  the  x-direction 
which  means  that  the  u-derivatives  there  are  larger  than  the 
v-derivatives.   It  follows  that  H  is  negative  and  so  &^   is 
negative  in  accordance  with  (15f).   The  same  reasoning  applies 
if  the  jet  length  is  too  small,  in  which  case  £^   is  positive. 
A  steady  state  is  reached  when  H  is  zero. 
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h.Ji.5'    Separation  Point 

The  separation  point  z   is  shifted  according  to  appro- 
priately applied  tangential  formulas,  with  the  constraint  that 
it  lies  on  the  line  x  =  0.   As  with  the  cusped  geometry  model, 
the  determination  of  the  proper  conformal  type  of  the  physical 
domain  is  an  important  part  of  the  problem.   For  a  rectangle 
In  the  w-plane  with  a  given  length-to-width  ratio,  we  expect 
to  achieve  a  conformal  map  onto  the  physical  domain  which 
satisfies  the  free  boundary  condition  only  if  the  separation 
point  is  properly  located  on  the  y-axis  relative  to  the  jet 
length.   In  the  cusped  geometry  model,  the  conformal  type  is 
established  by  the  normal  shift  term  only,  with  certain  con- 
straint conditions  along  the  walls  x  =  £,   y  =   £  as   discussed 
in  Section  4.2.3.   In  that  formulation,  the  apportionment  of 
points  between  fixed  and  free  boundaries  is  not  specified,  but 
is  determined  as  the  scheme  converges.   The  situation  is 
essentially  different  in  this  respect  for  the  vena  contracta 
model.   We  require  a  specific  point  on  the  rectangle,  namely 
the  corner  w  ,  to  map  into  the  separation  point  z^,  so  that  the 
number  of  points  which  map  onto  the  free  and  fixed  boundaries 
is  specified.   Since  z   is  shifted  sensibly  only  according  to 
tangential  terms  on  the  boundary,  we  see  that  the  normal  and 
the  tangential  shift  must  work  together  to  yield  the  proper 
conformal  type.   This  is  one  of  the  aspects  of  the  formulation 
of  the  vena  contracta  model  which  makes  the  issue  of  conver- 
gence considerably  more  delicate  than  it  is  for  the  cusped 
geometry  model. 
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The  formulas  which  we  use  to  shift  the  separation  point 
are  based  on  the  geometrical  interpretation  of  the  tangential 
shift  condition.   In  general,  a  steady  state  is  reached  when 
inner  normals  on  R  map  into  inner  normals  on  D  .   At  the 
separation  point,  where  the  right  angle  at  the  corner  w  is 
supposed  to  map  onto  a  l8o°  angle  at  z  ,  a  conformal  map  will 
carry  a  line  s  which  bisects  the  right  angle  in  R  onto  a  line 
s'  normal  to  the  l8o°  angle  in  D_.   Since  the  l8o°  angle  is 
oriented  in  the  vertical  direction,  this  normal  s'  is  a  hori- 
zontal line,  as  shown  in  Figure  4.6.   Therefore,  a  reasonable 
way  of  setting  up  the  tangential  shift  at  the  separation  point 
is  to  use  the  formulas 

X  =  0 
(13g) 


y^  =  a[(-y^  +  yv)//^]  • 


The  term  in  square  brackets  is  the  derivative  of  y  in  the  s- 
directlon  and  steady  state  is  achieved  in  this  formulation  when 
the  line  s'  becomes  horizontal,  characterizing  the  angle- 
doubling  property  of  the  conformal  map.   Another  way  of  looking 
at  these  formulas  is  to  observe  that  the  separation  point  can 
be  shifted  according  to  the  tangential  conditions  from  either 
the  fixed  boundary,  where  the  fact  that  x^  is  identically  zero 
yields  the  formulas 
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X  =  0 


or  from  the  free  boundary,  where.  If  we  assume  that  x  Is  zero 
which  It  Is  supposed  to  be  in  the  steady  state,  the  formulas 
are 

X  =  0 

yt  =  ^^^v  • 

Apart  from  the  factor  y?,  formulas  (13g)  represent  an  average 
of  these  tangential  conditions. 

We  argue  that  any  reasonable  method,  such  as  given  by 
(I5g),  for  treating  the  separation  point  which  results  in  an 
over-all  convergent  scheme  will  yield  a  finite  difference  solu- 
tion which  approximates  the  continuous  solution.   This  is  be- 
cause the  solution  to  the  continuous  problem  is  known  to  be 
unique,  and  in  particular,  there  is  only  one  conformal  mapping 
from  the  rectangle  onto  the  physical  domain  which  satisfies  the 
free  boundary  condition. 

We  obtain  an  additional  benefit  from  the  point  of  view 
of  stability  by  use  of  (l^g).   The  finite  difference  analog  of 
the  expression  -y  +  y  can  use  simple  forward  differences  to 
approximate  the  derivatives  instead  of  the  usual  three-point 
formulas  and  we  still  retain  second  order  accuracy  in  mesh 
size  h.   Specifically,  the  relation 
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-y^  +  y^  =  (y2- 2yQ+y-L)/h  +  o{h^) 


holds  In  the  steady  state  at  the  corner  w  ,  where  the  points 
are  labelled  as  shown  In  Figure  4.6.   This  follows  from  applica- 
tion of  a  Taylor  series  expansion  about  the  corner  and  use  of 
the  fact  that  the  solution  for  y  satisfies  Laplace's  equation 
at  the  corner  In  the  steady  state,  which  eliminates  the  first 
power  of  h  In  the  expansion.   The  use  of  three  points.  Instead 
of  five,  to  describe  the  motion  of  the  separation  point  as 
t  -►  00 ,  tends  to  Improve  the  stability  of  the  computation. 

4. 3.6.  Free  Boundary  Condition 

The  speed  q  at  any  point  along  the  free  boundary  r  Is 
given  by 

2    -2m, 2/,  2  ,  2, 
q  =  y   V^^v+^v^  • 

The  significance  of  A  In  the  normal  shift  term  of  equation  (13e) 
Is  that  Its  value  at  any  time  t  Is  an  estimate  of  the  constant 

speed  q  along  r,  and  so  In  what  follows  we  write  l/q  Instead 

2 
of  A  .   We  seek  the  steady  state  condition 


%/^   -1  =  0 


by  adjustment  of  the  free  boundary  normal  to  Itself  according 

to  equation  (l^e). 

o 
The  choice  for  l/q^  of  the  average 
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over  the  free  boundary,  which  Is  successful  In  the  cusped 
geometry  model,  yields  an  unstable  scheme  for  this  case,  as 
shown  by  the  following  heuristic  reasoning.   Suppose  the  free 
boundary  develops  a  corner  at  the  separation  point  as  the  tlme- 
dependent  problem  Is  being  solved.   Then,  under  the  assumption 
that  the  discrete  solution  at  this  stage  of  the  computation 
approximates  some  sort  of  flow,  the  speed  q  at  the  separation 
point  Is  high  relative  to  Its  average  value  q  over  the  free 


boundary,  and  we  have 


q^/q^  ^  1  <  0 


near  the  separation  point.   Since  the  normal  (x  ,y  )  Is  directed 
Into  the  physical  domain,  the  free  boundary  near  the  separation 
point  Is  shifted  outward,  tending  *to  Increase  the  sharpness  of 
the  corner.   The  same  argument  shows  that  a  bulge  In  the  free 
boundary  at  the  separation  point  will  also  be  accentuated,  hence 
Instability  of  the  free  boundary  at  the  separation  point  results 
from  the  averaging  formula. 

The  choice  suggested  by  the  above  analysis  of  Instability 
Is  to  make  l/q  correspond  to  the  speed  at  the  separation  point 
Itself.   Thus,  for  example.  If  a  sharp  corner  develops,  the 
separation  point  speed  q^  is  larger  than  the  speed  at  nearby  free 
boundary  points,  so  that  the  normal  shift  sends  these  points  into 
the  flow  region,  correcting  the  sharp  ctjrner.-  The  estimate  of 
q  must  be  obtained  as  a  limit  of  the  speed  along  the  free 
boundary  r  as  we  approach  the  separation  point  and  this  computa- 
tion Is  carried  out  numerically  by  fitting  a  linear  function  to 
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values  of  l/q  nearby  and  extrapolating  to  the  separation 
point.   The  slope  of  this  linear  function  must  go  to  zero  as 
t  -+  00  .   When  the  slope  is  different  from  zero,  it  essentially 
provides  a  stabilizing  correction  to  the  averaging  formula. 

We  must  retain  in  the  normal  shift  formula  the  provision 
for  making  the  conjugate  function  x^+y^  -  x^  -  y^  equal  to  zero 
in  the  steady  state.   The  general  theory  shows  that  this  func- 
tion is  equal  to  a  constant  H  when  the  Dirichlet-Douglas  func- 
tional is  minimized,  so  that  if  H  =  0  at  even  one  point  on  the 
boundary,  for  example  the  separation  point,  then  it  is  zero 
throughout  R.   To  achieve  the  condition  H  =  0,  we  must  incorpo- 
rate u-derivatives,  which  are  tangential  derivatives  on  r,  into 

the  normal  shift  formula.   We  recall  that  in  the  cusped  geometry 

2    2 
model  this  is  achieved  by  use  of  an  average  of  x^  + y^  over  r  so 

that  as  t  -»•  00  the  condition  H  =  0  is  obtained  numerically,  not 

Just  at  one  point,  but  over  the  whole  free  boundary.  Although 

2    2 
the  theory  suggests  that  we  could  use  the  value  of  x^+y^  at 

Just  one  point  on  the  free  boundary  in  the  normal  shift  formula, 

such  a  formulation  is  not  expected  to  produce  a  stable  scheme. 

The  situation  is  analogous  to  solving  an  elliptic  equation  in 

which  we  specify  a  directional  derivative  with  a  non-zero 

tangential  component,  that  is,  a  u-derivative  along  the  boundary, 

and  such  problems  are  difficult  to  treat.   The  use  of  averaging, 

equivalent  to  integration  along  the  boundary,  is  one  way  out 

of  the  difficulty.   For  the  vena  contracta  model,  we  have  shown 

that  this  does  not  work.   We  are  forced  to  use  a  value  of  l/q^ 
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based  on  a  single  point,  and  u-derlvatlves  in  the  extrapolation 
cause  instability.   Therefore  we  must  use  v-derlvatives,  which 
are  normal  derivatives  on  the  free  boundary,  to  estimate  l/q 
at  the  separation  point,  and  we  are  still  faced  with  the  task 
of  building  into  the  normal  shift  formula  a  mechanism  which 
yields  H  =  0  in  the  steady  state.   We  do  this  by  numerically 
Imposing  the  requirement  that  the  speed  be  continuous  at  the 
separation  point.   Continuity  of  the  speed  is  implicit  in  the 
continuous  formulation  of  the  problem,  but  not  in  the  discrete 
formulation,  where  it  must  ultimately  be  allied  with  the 
conformallty  of  the  mapping. 

Our  requirements  lead  us  to  the  specification  for  l/q 
as  the  weighted  average  given  by 

(16)  l/q^  =  (l-e)/q^  +  e/q^ 

where  e  is  the  weight  parameter,  chosen  experimentally  to 
produce  convergence,  and  q-,  and  q^  are  the  speeds  extrapolated 
from  the  free  boundary  and  the  fixed  boundary,  respectively, 
given  by  the  limits 

1/q^  =  lim  [y^"^(x2+y2)/^2^ 
v=0 

1/q^  =  lim  [y2"^(x^  +  y2)/^2^  . 
u=a 
v|0 

These  limits  are  evaluated  numerically  by  extrapolation  of  a 
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linear  function  which  fits  the  data  In  a  least  squares  sense 
at  a  suitable  number  of  points.   On  the  free  boundary,  the 
truncation  error  Is  not  affected  by  the  number  n^  of  points 
used.   In  steady  state  the  slope  of  the  linear  function  Is 
zero  because  the  speed  Is  constant  all  along  the  free  boundary. 
The  number  n  Is  therefore  chosen  solely  on  the  basis  of 
stability.   On  the  fixed  boundary,  the  speed  Is  not  constant 
so  that  the  number  n^  of  points  used  In  the  linear  extrapola- 
tion should  be  taken  as  small  as  possible  consistent  with 
stability.   Of  course  we  need  a  minimum  of  2  points  In  order 
to  define  a  slope. 

The  conditions  that  the  speed  Is  continuous  and  that  H 
Is  zero  In  steady  state  are  achieved  In  the  following  way. 
Prom  the  steady  state  normal  shift  term  In  (13e)  we  have 

along  r,  so  that  the  free  boundary  condition  Is  satisfied.   As 
a  consequence,  we  must  also  have 

l/q^=l/q^ 

Since  l/q^  Is  an  extrapolation  of  the  (constant)  values  of 
l/q^  along  the  free  boundary.  Therefore,  by  Inspection  of 
equation  (l6)  we  also  have 

1/q^  =  1/q^  , 
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so  that  the  speeds  q, ,  qp  and  q  are  all  equal.   To  show  that 
H  Is  zero  we  argue  that  since  ip     and  f     both  go  to  zero,  and 

the  limits  l/q-,  and  l/q^  are  finite  at  the  separation  point, 

2   2      2   2 

then  it  must  follow  that  the  quantities  x  +y  and  x  +y  , 

^  V  "'v      u   u' 

which  are  the  numerators  in  these  limits,  also  both  go  to  zero 
there.   In  this  way  we  obtain 


H  =  0 

at  the  separation  point,  and  by  the  general  theory,  the  mapping 
is  conformal.   We  remark  that  there  is  an  analytic  reason  why 
the  weight  factor  e  should  have  to  be  chosen  small.   If  there 
exist  solutions  to  the  differential  equations  which  satisfy  all 
steady  state  conditions  with  the  exception  that  the  speed  is 
not  necessarily  continuous  at  the  separation  point,  then  these 
solutions  have  the  property  that  both  x  and  y  are  zero  at  the 
separation  point  because  the  limit  of  speed  coming  from  the 
free  boundary  is  finite.   Also,  by  the  way  we  treat  the  tangen- 
tial shift  of  separation  point,  the  derivative  y  is  zero. 
There  is  no  control,  though,  on  x  ,  which  can  be  different 
from  zero.   The  result  is  that  l/qo^  the  limit  from  the 
fixed  boundary  side  of  the  separation  point,  is  infinite,  due 
to  the  non-zero  x  .   A  small  e  is  necessary  to  prevent  in- 
stability  caused  by  the  development  of  large  values  of  l/qp  as 
the  scheme  converges. 
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4.3.7.  Rlabouchinsky  Model 

The  application  of  the  numerical  method  to  the 
Rlabouchinsky  model  is  the  same  as  it  is  for  the  vena  contracta, 
and  so  we  shall  omit  the  details.   Here  also,  by  mapping  a 
corner  of  the  rectangle  into  the  separation  point,  we  uniformize 
the  singularity  in  the  stream  function,  and  make  the  over-all 
scheme  accurate  to  second  order  in  mesh  width  h.   The  following 
observations  are  noted: 

1)  The  tangential  shift  of  the  separation  point,  calculation 
of  the  speed  as  a  weighted  average  of  the  limits  from  the  free 
and  fixed  boundaries,  the  introduction  of  an  interface  near 
infinity  and  subtraction  of  the  pole  in  the  desired  conformal 
mapping,  and  the  associated  horizontal  shift  in  the  vertical 
line  of  symmetry  of  the  flow,  which  are  techniques  we  have 
described  in  detail  for  the  vena  contracta  model,  are  carried 
over  and  applied  in  exactly  the  same  way  for  the  Rlabouchinsky 
model. 

2)  The  stream  function  can  be  written  in  the  form 

f  =   Uy^"^"'/(l+m)  +  ^j^ 


where  the  regular  part  t^„  -♦  0  as  |z|  — »•  oo.   In  the  triangular 
subdomain  R,  defined  by  the  diagonal  interface  we  solve  numeri- 
cally for  Tp^   in  the  plane  case  (m  =  0),  and  for  f^   in  the 

K  K 


axlally  symmetric  case  (m  =  1),  where 


^R=  1^1% 
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Here  the  Kelvin  inversion  (Section  4.1.5)  takes  the  slightly 
different  form 

to  account  for  the  -g-shift, 

3)  The  boundary  correspondences  for  the  conformal  map  are 
in  accordance  with  those  given  in  Figure  2.3a  and  2.3c,  for 
both  the  plane  and  axially  symmetric  models.   The  initial 
mapping  is  taken  in  either  of  two  ways.   In  the  plane  model  it 
is  the  exact  solution  to  the  plane  flow  model  given  by  equation 
(2.6),  for  a  specific  value  of  k.   It  is  possible  to  choose  k 
so  that  the  dimensions  of  the  rectangle  K  and  K  have  a 
rational  ratio,  allowing  for  equal  mesh  spacing.   In  the  axially 
symmetric  model,  the  plane  solution  is  used  as  the  initial 
mapping  for  crude  mesh  widths,  and  for  finer  meshes  an  inter- 
polated solution  is  used. 

4.3.8.  Remarks 

We  view  the  following  features  as  the  three  delicate 
aspects  of  the  numerical  method: 

1)  The  equations  (13g)  for  shifting  the  separation  point. 

2)  The  specification  of  l/q_  by  equation  (16)  in  the  normal 
shift  term. 

3)  The  treatment  of  the  neighborhood  of  infinity  and  the 
associated  ^-shift  (13f)  of  the  vertical  line  of  symmetry. 

The  formulas  and  procedures  used  were  arrived  at  for  the 
most  part  experimentally.   Taken  individually,  there  are  many 
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reasonable  ways  in  which  to  deal  with  each  of  these  aspects. 
However,  in  practice  we  found  that  it  was  a  more  difficult 
task  to  make  them  all  work  together.   It  is  not  yet  entirely 
clear  why  some  ideas  failed  and  the  ones  we  have  given  were 
successful.   Moreover,  it  is  quite  probable  that  still  other 
approaches  can  greatly  improve  the  convergence  properties 
compared  to  what  we  have  achieved  to  date. 

We  include  a  final  remark  on  the  heuristic  explanations 
of  several  aspects  of  the  scheme.  We  use  formulas  which  make 
the  quantity 

2,2    22 

H  =  x  +y   -X  -y 
V   "^v    u  ''u 

equal  to  zero  both  at  the  separation  point  and  in  the  ^-shift 
formula.   This  is  really  contradictory,  since  the  theory  says 
that  when  we  have 

X  X  +  y  y  =0 
u  V  •^u-'v 

on  SR,  then  the  harmonic  conjugate  H  is  identically  constant 
throughout  R.   A  more  sophisticated  theory  takes  into  account 
the  fact  that  the  analytic  function  z(w)  has  a  pole  at  one  of 
the  corners  of  R,  the  treatment  of  which  necessitated  the 
addition  of  the  ^-shift  (l^f)  in  the  first  place.   The  harmonic 
function  ^vi-^y  +  y^nYy  also  has  a  pole  at  this  corner,  which  means 
that  it  is  possible  that  it  can  be  zero  everywhere  else  on  the 
boundary  of  R,  without  being  zero  in  the  interior.   In  con- 
sequence, the  normal  derivative  of  the  harmonic  conjugate  H  is 
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zero  on  c)R,  but  H  is  not  necessarily  constant  in  R.   This  is 
analogous  to  the  example  of  the  function  f{z)  =  z,  whose 
Imaginary  part  y  is  zero  on  the  x-axls,  and  positive  in  the 
upper  half -plane,  and  whose  real  part  has  a  vanishing  normal 
derivative  on  the  x-axls,  and  a  nodal  line  (x  =  0)  extending 
to  Infinity.   In  our  numerical  scheme,  we  impose  the  condition 
H  =  0  at  the  separation  point  as  described  in  Section  h.J).G. 
This  means  that  we  obtain  a  nodal  line  (H  =  0)  which  emanates 
from  the  separation  point  and  runs  into  R.   We  reason  that 
this  nodal  line  must  not  be  allowed  to  wander  off  to  the  corner 
at  which  the  pole  is  located,  for  by  analogy  with  the  above 
example  of  f(z)  =  z,  such  behavior  means  a  non-unique  solution 
(cf.  the  Phragman-Lindelof  principle).   We  therefore  use  the 
^-shift  formula  to  control  the  path  of  the  nodal  line  so  that 
it  terminates  on  some  other  part  of  SR,  thereby  enclosing  some 
finite  subdomain  of  R,   If  this  is  the  case,  then  we  can  use 
well-known  theorems  to  conclude  that  H  =  0  in  this  subdomain. 
Prom  this  it  follows,  by  analytic  continuation,  that  the 
analytic  function  for  which  H  and  x^x^  +  y^y^  are  the  real  and 
imaginary  parts,  is  identically  zero  throughout  R,  and  we 
obtain  the  Cauchy-Riemann  equations  for  z(w). 
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5.  RESULTS 

The  main  result  is  that  the  numerical  method,  when 
applied  to  the  models  we  have  studied.  Is  convergent  for 
appropriately  chosen  factors.   Results  for  the  cusped  geometry 
model  Indicate  the  expected  high  degree  of  stability  except  in 
geometric  configurations  for  which  the  effects  of  truncation 
error  at  the  separation  point  become  pronounced.   For  the  hydro- 
dynamic  models,  both  plane  and  axlally  symmetric  calculations 
are  performed  for  mesh  widths  sufficiently  small  to  demonstrate 
that  the  numerical  solution  converges  to  the  continuous  solu- 
tion.  Computations  show  that  the  contraction  coefficient  and 
drag  coefficient  for  the  plane  models  agree  to  at  least  four 
significant  figures  with  the  exact  values,  and  new  estimates 
are  obtained  for  these  quantities  In  the  axlally  symmetric  case. 
Calculations  are  also  described  In  which  the  method  Is  applied 
as  a  technique  for  conformal  mapping  only. 

5.1.  Cusped  Geometry  Model 

The  geometric  configuration  for  this  model  Is  controlled 
by  specification  of  the  wall  constraint  I   and  the  radius  r  of 
the  coil.    In  practice,  r  was  fixed  by  the  height  b  of  the 
rectangle  according  to  the  relation 

-b 
r  =  e 

which  scales  the  physical  domain  so  that  ^  must  be  selected  in 
the  range  v  <   I   <   \   for  non-trivial  solutions. 
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The  values  of  the  convergence  factors  a,    p,  and  7  were 
determined  experimentally  to  produce  rapid  convergence  con- 
sistent with  the  stability  of  the  computation.   The  value  of 
a  must  be  sufficiently  large  to  provide  damping  of  the  interior 
equations.   Selection  of  a  is  equivalent  to  setting  the  value 
of  the  relaxation  factor  in  these  equations.   The  factors  P 
and  7  control  the  size  of  the  motion  of  the  boundary  points. 
In  accordance  with  the  theory  of  the  method  of  steepest  descent, 
P  and  7  must  be  sufficiently  small  to  prevent  overshooting 
the  minimum  point  of  the  multi-dimensional  surface  given  by 
the  numerical  formulation. 

Table  I  and  Figure  5«1  summarize  the  results  of  a 

specific  calculation.   Our  purpose  in  conducting  this  sequence 

of  experiments  is  to  examine  the  convergence  properties  of  the 

iterative  solution  to  a  steady  state  answer,  where  steady  state 

_o 
is  defined  as  the  reduction  of  all  residuals  to  less  than  10 

We  also  want  to  observe  the  behavior  of  the  discrete  solution 

as  the  mesh  is  successively  refined.   We  note  the  following 

remarks. 

1.  The  number  of  iterations  required  to  reach  steady  state 
approximately  doubles  with  each  refinement  of  the  mesh  by  a 
factor  of  two.   This  confirms  the  linear  theory  described  in 
Section  4.2.3  which  asserts  that  the  number  of  iterations  grows 
like  0(l/h). 

2.  About  19  micro-seconds  are  required  to  perform  an 
iteration  in  both  x  and  y  at  a  single  interior  mesh  point  on 
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the  CDC  6600  computer.   At  this  speed  the  problem  for  the  crude 
mesh  h  required  about  9  seconds  to  converge,  and  the  problem 
for  the  most  refined  mesh  h/8  required  about  56  minutes. 

3.  Figure  5.1  shows  the  computed  free  boundary  for  this 
geometric  configuration  using  the  data  for  the  hA  mesh.   In 
this  solution  the  walls  r^  consist  of  19  points  each  and  the 
free  boundary  r  is  described  by  75  points.   The  tabulation  of 
the  separation  points  S^  and  Sg  and  the  midpoint  M  in  Table  I 
for  the  various  mesh  widths  indicates  the  expected  behavior 
due  to  truncation  error.   The  midpoint  location  is  stabilized 
to  three  digits  at  (.719,-719)  but  the  location  of  the  separa- 
tion points  is  known  only  to  be  in  the  range  .21-. 22.   This 

is  because  no  provision  was  made  to  deal  with  the  known 
singular  behavior  at  these  points. 

4.  The  cusped  geometry  model  provides  an  example  with 
a  curved  fixed  boundary,  namely  the  circular  coil  O,  in  which 
second  order  normal  corrections  are  made  to  keep  tangentially 
shifted  points  on  the  boundary,  as  described  in  Section  4.2.5. 
The  procedure  was  quite  satisfactory  in  this  model  and  we 
expect  that  curved  boundaries  in  general  can  be  treated  in  the 

same  way. 

5.  The  initial  mapping  used  for  most  runs  was  the  trans- 
formation described  in  Section  4.2.2.   In  order  to  verify  that 
the  convergence  is  not  one-sided,  runs  were  also  made  in  which 
the  initial  mapping  specified  a  free  boundary  T    close  to  the 
corner  &+±£    so  that  the  motion  toward  steady  state  was  inward. 
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away  from  the  corner.   We  were  able  to  obtain  convergence  for 
a  wide  variety  of  such  Initial  conditions,  demonstrating  the 
high  degree  of  stability  possessed  by  the  scheme. 

Experiments  were  also  made  In  which  the  value  of  the  wall 
constraint  H   was  reduced,  requiring  the  free  boundary  r  to  move 
closer  to  the  cover  i +  li  In  order  to  satisfy  the  condition  of 
conf ormallty.   This  effect  reduces  the  number  of  points  on  the 
base  of  the  rectangle  which  map  onto  r.   The  results  of  these 
calculations  are  qualitatively  similar  to  what  we  have  described 
In  Table  I.   As  r  Is  reduced,  however,  the  truncation  error  In 
the  vicinity  of  the  separation  builds  up  rapidly  due  to  the 
Increased  stretching  of  the  physical  domain.   There  Is  a  fairly 
sharp  threshold  In  S, ,    below  which  the  free  boundary  appears  to 
"roll  up"  toward  the  corner  and  Is  numerically  described  by  just 
one  point,  the  midpoint  M,  with  all  the  other  points  pasted  to 
the  walls  T^,      Refinement  of  the  mesh  Is  necessary  to  retain 
accuracy  of  the  numerical  solution.   The  same  effect  Is  observed 
when  calculations  are  made  with  smaller  values  of  the  radius  r. 
This  behavior  exemplifies  the  Importance  of  treating  the  region 
near  the  separation  point  properly  If  a  high  degree  of  accuracy 
is  to  be  attained. 

Calculations  were  also  made  in  which  the  wall  constraint 
was  removed,  that  is,  H   >   1.      In  this  limiting  case  the  separa- 
tion points  coincide  with  the  corners  of  the  outer  boundary  and 
are  not  singular  points.   In  effect,  this  model  is  a  variant  of 
the  trivial  rectangle  mapping  described  in  Section  3-5' 
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Experiments  verify  that  the  convergence  of  the  numerical  solu- 
tion to  the  continuous  trivial  solution  z  =  e~   is  0(h  ). 

5.2.  Hydrodynamic  Models 

In  this  section  we  describe  the  results  of  the  calcula- 
tions of  the  vena  contracta  and  Riabouchinsky  flows.   These 
are  models  of  the  type  for  which  the  method  was  designed  to  be 
accurate  to  O(h^)  by  mapping  a  corner  of  R  into  the  separation 
point.   Our  general  results  are  as  follows: 

1.  The  scheme  is  convergent  for  a  wide  variety  of 
geometrical  configurations  provided  that  the  length  to  width 
ratio  r  of  the  parameter  rectangle  R  is  not  too  large.   By 
taking  the  mesh  width  h  sufficiently  small  we  verified  the 
convergence  to  be  0(h  ). 

2.  The  method  works  equally  well  for  axlally  symmetric 
and  plane  geometries. 

3.  Drag  and  contraction  coefficients  for  the  plane  model 
agree  with  their  exact  values  to  almost  four  and  almost  five 
significant  figures,  respectively.   For  the  axially  symmetric 
models,  new  estimates  of  these  quantities  are  obtained  which 
are  judged  to  be  as  accurate  as  the  results  of  the  calculations 
for  the  plane  models. 

The  details  are  given  in  the  following  sections. 

5.2.1.  General  Results 

Calculations  of  the  hydrodynamics  models  were  made  for 
various  mesh  sizes  h  and  length-to-height  ratios  r  of  the 
parameter  rectangle.   In  Figure  5-2  we  give  the  results  of 
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particular  calculations  in  terms  of  graphs  of  the  profiles  of 
the  free  boundary  for  the  plane  and  axially  symmetric  jet.   In 
Figure  5.3,  the  same  kinds  of  results  are  shown  for  the 
Riabouchinsky  model.   We  observe  the  extremely  fine  resolution 
near  the  separation  point,  as  shown  in  the  enlarged  inserts, 
where  the  actual  data  points  are  plotted.   The  high  density  of 
mesh  points  in  this  region  is  the  consequence,  numerically,  of 
having  mapped  a  corner  of  R  into  the  separation  point.   For  both 
models,  the  exact  plane  solution  overlays  the  numerically  com- 
puted solution  for  the  resolution  shown. 

In  practice,  the  various  convergence  factors  which  occur 
in  the  formulation  of  these  problems  were  determined  by  trial 
and  error  to  produce  convergent  schemes.   These  factors  are 
listed  below: 

1)  a  -  dissipation  factor  for  the  interior  equations,  which 
controls  the  damping  (relaxation  factor)  in  equation  (4.1), 

2)  P,  7  -  convergence  factors  for  the  tangential  and  normal 
shifts  at  the  boundary. 

3)  6,  Re  -  convergence  factor  for  the  i-shift  of  the  line 

of  symmetry  and  the  subdomain  Rg  over  which  the  average 

in  (4.13f)  is  computed. 

i\)      n  ,  n'  -  number  of  points  used  in  the  extrapolation  of  the 
e   e 

speeds  from  the  free  and  fixed  boundaries,  respectively. 
5)   i   -  the  number  of  points  along  the  vertical  side  of  the 
rectangle  above  which  the  triangular  subdomain  R^  is 
defined.   This  parameter  locates  the  intersection  of  the 
diagonal  interface  with  the  vertical  side  of  the  rectangle, 
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6)   e  -  weight  factor  In  the  computation  of  the  weighted 
average  of  speeds  at  the  separation  point. 

In  Table  II  we  list  the  appropriate  data  for  the  specific  cases 
we  have  given.   Calculations  were  also  made  for  various  other 
values  of  r  and  h  in  order  to  obtain  data  for  the  contraction 
coefficients  and  drag  coefficients.   The  following  observations 
are  pertinent  for  these  runs: 

1)  Verification  that  the  numerical  solution  converges  with  h 
to  the  continuous  solution  is  given  in  Table  III.   Data  is 
shown  for  a  particular  axially  symmetric  vena  contracta  and 
plane  Riabouohinsky  model.   The  quantities  tabulated  are  the 
location  of  the  line  of  symmetry  &   and  Its  u-derlvatlve,  which 
are  typical  data  points.   Calculations  show  that  the  convergence 
is  0(h  )  for  sufficiently  small  h. 

2)  The  success  of  the  scheme  was  fairly  sensitive  to  the 
choice  of  an  initial  mapping.   In  practice,  for  a  calculation 
with  a  refined  mesh,  we  found  It  best  with  the  axially  symmetric 
models  to  use  as  an  initial  mapping  the  solution  previously 
calculated  for  a  cruder  mesh  Interpolated  at  the  new  mesh  points. 
In  some  cases  it  was  found  expedient  to  heavily  damp  the  time- 
dependent  equations  during  the  first  few  hundred  cycles  of 
iteration  in  order  to  smooth  out  the  data.   The  number  of  itera- 
tions required  for  convergence  obeys  the  linear  theory  of 
Section  4.1.5  for  the  axially  symmetric  vena  contracta  as  shown 
by  the  data  of  Table  III.   On  the  other  hand.  Table  III  indicates 
a  faster  rate  of  convergence  for  the  plane  Riabouchlnsky  model. 
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This  is  because  as  the  mesh  is  refined,  the  initial  mapping 
provides  better  initial  data  so  that  the  numerical  solution 
starts  closer  to  the  final  answer. 

5)   Experiments  show  that  to  iterate  1000  mesh  points  for 
1000  cycles  takes  about  8o  seconds  for  the  axially  symmetric 
case  and  about  35  seconds  for  the  plane  case  on  the  CDC  6600 
computer.   In  practice,  some  economy  was  achieved  due  to  the 
fact  that  the  convergence  of  the  f   equation  was  quite  rapid. 
This  permitted  iteration  of  the  f   equation  only  once  for  every 
several  cycles  of  iteration  of  the  x  and  y  equations.   For 
the  vena  contracta  model  we  were  able  to  use  a  ratio  of  1  to 
10  in  this  scheme.   For  the  Riabouchinsky  model,  where  the 
i-shift  affects  the  i{/   equation  more  directly,  ratios  of  1  to 
5  were  possible  for  the  lower  valties  of  r,  but  we  had  to  use 
a  1  to  1  ratio  for  the  largest  value  r  =  25/12.   The  actual 
computing  time  to  reduce  residuals  below  10"  was  about  40 
minutes  for  the  plane  vena  contracta  model  and  about  65  minutes 
for  the  axially  symmetric  vena  contracta  model  of  Table  II. 

4)  Values  of  the  weight  factor  e  in  the  range  .02-.04  for 
the  Riabouchinsky  model  and  about  .1  for  the  vena  contracta 
model  were  found  by  experiment  to  produce  convergent  schemes. 
This  confirms  the  analysis  given  in  Section  4.3.6  that  small 
values  of  e  are  required  in  the  formula  for  l/q^. 

5)  The  location  of  the  interface,  governed  by  the  parameter 
£    ,    was  determined  experimentally,  but  also  with  the  motivation 
that  the  derivatives  along  it  should  be  reasonably  close  to 
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unity,  so  that  the  transmission  of  truncation  error  between 
the  regions  is  minimized.   In  practice,  we  found  that  a  variety 
of  interface  locations  were  acceptable  to  produce  convergence. 

6)  Values  of  a,    p,  7,  and  5  were  found  experimentally. 

7)  For  the  vena  contracta  model,  the  region  Re,  over  which 
the  average  in  equation  (4,13f )  is  computed,  was  taken  to  be 
the  entire  rectangle  excluding  the  triangular  region  R, .   In 
the  Riabouchinsky  model,  experiments  with  large  r  have  shown 
that  a  different  choice  of  R^  is  necessary  to  produce  conver- 
gence.  We  choose  Rg  to  be  a  fixed  segment  MA'  along  the  base 
MA  of  the  rectangle  (Figure  2.3c).   Furthermore,  we  normalize 

0000 

(4.13f)  by  dividing  each  term  in  the  average  by  x  +y  +  x  +y 
so  that  all  terms  have  equal  weight.   In  view  of  the  nodal  line 
theory  developed  in  Section  4.3.8,  we  argue  that  this  recipe 
exerts  a  stronger  influence  in  making  the  nodal  line  exit  from 
the  interior  of  R  at  a  point  along  the  base  MA,  whereas  a 
choice  of  R^  as  in  the  vena  contracta  model  permits  the  nodal 
line  to  wander  off  to  the  corner  where  the  pole  is  located. 
Further  experiments  are  in  progress. 

5.2.2.  Contraction  Coefficient 

We  used  our  numerical  method  to  calculate  the  contraction 
coefficient  C^  for  the  infinite  jet  (r -►  00).   The  known  value 
in  the  plane  case  (m  =  0),  is  given  by 


^o  =il?=  • 611015  , 
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and  can  be  used  to  check  the  computation.   In  the  axlally 
symmetric  case  (m  =  1),  values  of  C^  In  the  range  .57-  .6l 
have  resulted  from  various  other  techniques  of  estimation 

[5,6,13,1^]. 

By  using  as  long  a  rectangle  and  as  small  a  mesh  width 
as  possible  we  find  that  the  computed  value  of  C  is 


Cq  =  .61129 


which  shows  agreement  with  the  exact  answer  to  three  digits. 
However,  we  can  improve  our  estimate  by  making  use  of  the 
expected  asymptotic  behavior  of  C  as  h  -»■  0  and  r  -»■  oo  and 
performing  suitable  extrapolation  of  the  data.   In  Section  2.5 
we  obtained  a  formula  for  this  behavior,  given  by 


-A  r 

(1)  C„(r)  -  C^+A  e  ^ 

^  '  m      m   m 


as  r  -►  00  ,  where  the  constant  A  has  the  values 

'  m 


■K      =  5  ^  1.6 
o    2 


A^  =  2.4   . 


We  denote  by  C  (r,h)  the  numerically  computed  contraction 
coefficient  for  a  given  r  and  h.   Since  the  numerical  scheme 

is  0(h  )  this  quantity  should  converge  for  fixed  r  with  an 

2 
error  term  M  h  .   The  number  M  differs  from  a  constant  by 
m  m 

terms  which  go  to  zero  exponentially  with  r.   Likewise,  for 


90 


-A  r 
fixed  h,  C  {r,h)  converges  with  an  error  term  A  e  '^  ,  where 

A  and  A  differ  from  constants  by  terms  which  are  0(h  ). 
mm 

Therefore,  we  expect  the  asymptotic  behavior  of  C  (r,h)  to  be 
given  by 

(2)  C^{r,h)  -C^+M^+A^e  ^"^ 


as  r  -*■  00  ,  h  — ►  0,  where  A  ,  A  and  M  are  assumed  constant. 

'  '  m'   m      m 

To  carry  out  our  estimate  we  made  a  series  of  computations  of 

C  (r.h)  which  are  summarized  in  Table  IV  and  Table  V.   The 
m 

columns  of  these  tables  indicate  the  ratio  r.   The  rows  indicate 

mesh  refinements  in  fractions  of  the  basic  width  h  =  l/l2.   We 

remark  that  the  case  r  =  6,  3/l6  h  in  Table  IV  corresponds  to 

25,025  mesh  points,  which  is  an  upper  limit  in  terms  of  core 

storage  capacity  on  the  CDC  660O  computer. 

According  to  equation  (2)  the  coefficients  A_  and  A 
°     ^       ^  '  mm 

should  be  independent  of  which  row  we  use  in  our  tables.   Our 
calculations  confirm  this  and  we  obtain 

A^  =  1.9  ,    A^  =  1.6  , 


A^  =  2.3  ,    A^  =  2.4  . 


The  agreement  of  A  with  our  estimates  is  correct  to  two 
•^  m 

significant  digits. 

We  remark  that  the  axial ly  symmetric  jet  is  expected  to 
converge  with  r  much  faster  because  of  the  relative  sizes  of 
A^  and  A  .   Our  calculations  confirm  this  behavior.   Using 
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equation  (2),  the  correction  for  a  jet  corresponding  to  r  =  4.5 
is  in  the  fifth  decimal  place.   This  is  why  we  only  had  to  take 
r  =  4.5  for  the  axially  symmetric  model,  whereas  for  the  plane 
jet  a  calculation  for  r  =  6  was  required. 

For  the  extrapolation  with  h,  we  observe  that  appreciable 

2 

refinement  of  the  mesh  is  necessary  for  the  term  Mh  to  dominate 

the  error  as  h  — >•  0.   In  Table  V,  calculations  show  that  the 
convergence  is  0(h  )  for  the  h/8  mesh  (r  =  1.5)  and  for  the 
3/16  h  mesh  (r  =  3.75).   For  the  jet  lengths  corresponding  to 
r  =  3  and  r  =  ^-5,   where  the  most  refined  mesh  was  only  h/4,  we 
find  that  there  are  contributions  to  the  error  from  higher 
order  terms.   The  same  effect  is  present  in  the  plane  model, 
where  we  had  to  go  to  a  3/16  h  mesh  (r  =  6)  to  verify  that  the 
convergence  is  0(h  ). 

Calculations  using  the  extrapolation  formula  (2)  give  the 
final  estimates 

Cq  =  .61100 
c^  =  .59135 


The  value  of  C  is  in  error  by  only  2  units  in  the  fifth  decimal 
place.  Close  analysis  of  the  data  leads  us  to  conclude  that  the 
error  tolerance  in  the  value  of  C-,  is  ±.00004, 

5.2.3.  Drag  Coefficient 

The  drag  coefficient  C^(a)  was  computed  from  the  finite 
difference  solution  of  the  Riabouchinsky  model  by  the  formula 
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y 


Cj^Ca)  =  (1+a) 


K' 


0 


which  Is  easily  derived  from  the  basic  definition  of  C^.   In 
this  formula  ct  is  the  cavitation  parameter,  q  is  the  computed 
speed  along  the  free  boundary,  y^  is  the  height  of  the  plate, 
m  is  0  or  1  for  the  plane  or  axially  symmetric  case,  respec- 
tively, and  the  integration  is  carried  out  along  the  vertical 
side  of  the  rectangle  which  maps  onto  the  obstacle  according  to 
Simpson's  Rule. 

As  with  the  vena  contracta  model,  we  calculated  a(r,h) 
and  Cj^(r,h)  for  various  values  of  the  rectangle  length-to-width 
ratio  r  =  K/K'  and  for  the  three  levels  of  mesh  refinement  h, 
h/2,  h/4,  where  h  =  K'/l2  (i.e.  12  intervals  in  the  vertical 
direction).   The  results  of  the  computations  are  shown  in 
Tables  VI  and  VII  for  the  plane  and  axially  symmetric  cases, 
respectively. 

For  the  plane  case,  extrapolation  with  h  yields  the  data 
shown  in  the  next  to  last  row  of  Table  VI.   By  comparison  with 
the  last  row  of  the  Table,  which  shows  the  exact  values  of  o 
and  Cp/l+a  for  the  given  values  of  r,  we  see  that  our  computa- 
tions are  accurate  to  about  4   significant  figures.   The 
maximum  value  13/6  for  r,  corresponding  to  a( exact)  =  .305355 
represents  an  upper  limit  on  r  in  terms  of  convergence,  for  we 
found  that  for  larger  values  the  scheme  diverged.   This  diver- 
gence, which  of  course  occurs  also  in  the  axially  symmetric  case, 
is  the  inevitable  consequence  of  the  rapidly  increasing  ratio 
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of  cavity  length  to  obstacle  length  caused  by  using  the 
Riabouchlnsky  model  to  approximate  Helmholtz  flow.   The  free 
streamlines  for  Helmholtz  flow  actually  have  a  parabolic  shape 
at  large  distance  from  the  obstacle. 

For  the  axially  symmetric  case,  extrapolation  with  h  yields 
the  data  shown  in  the  last  row  of  Table  VII.   We  estimate  the 
accuracy  to  be  comparable  to  that  obtained  in  the  plane  case. 
The  value  r  =  25/12  is  an  upper  limit  for  convergence  of  the 
scheme.   The  results  are  plotted  as  a  graph  in  Figure  5.4.   We 
observe  that  the  graph  is  remarkably  linear  in  the  range 
.08  _<  a  ^  .6  and  that  for  smaller  a   there  is  an  indication  of 
flattening  of  the  curve.   This  qualitatively  supports  the  theory 
[5]  for  the  axially  symmetric  model  which  states  that  the 
asymptotic  variation  of  C^/l+a  is  o(a).   Furthermore,  extrapola- 
tion to  a  =  0  yields,  to  5  significant  digits,  the  result 

Cj^{0)  =  .828  , 

which  is  in  good  agreement  with  the  estimate  stated  in  [5]^ 

c^(o)  =  .827  . 

Finally  we  remark  that  our  results  are  reasonably  placed  on 
the  graph  of  experimental  data  available  in  the  physically 
significant  range  .1  ^  o  ^  .J>* 

5.3.  Conformal  Mapping 

As  stated  in  Section  3.2,  the  method  of  steepest  descent 
applied  to  the  minimization  of  the  Dlrlchlet-Douglas  functional 
provides  a  numerical  technique  for  obtaining  conformal  mappings. 
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In  this  case,  a  3-point  condition  is  specified  and  only 
tangential  shifts  are  applied.   In  the  course  of  our  investiga- 
tion of  the  more  difficult  free  boundary  problem,  we  conducted 
experiments  using  the  conformal  technique  alone.   The  results 
are  fragmentary,  but  they  indicate  that  the  method  can  be 
applied  quite  successfully. 

We  compute  the  analytic  function  which  maps  the  rectangle 
R  conformally  onto  the  physical  domain  D  as  shown  in  Figure 

5.5. 

The  following  remarks  are  pertinent: 

1)  The  rectangle  R  is  fixed  as  shown  in  the  w-plane.   The 
image  domain  D  consists  of  the  part  of  the  first  quadrant 

bounded  by  the  segment  x  =  b,  Oj^yj^a,  an  arc  of  the  circle 

2       2    2 
X  +  (y-b)   =  a  ,  and  appropriate  segments  of  the  coordinate 

axes. 

2)  The  3-point  condition  is  given  for  the  corners  of  the 
rectangle  labelled  A,  C,  and  D.   In  particular  the  point  D  is 
mapped  into  00  in  the  z -plane.   The  corner  B  is  free  to  move 
along  the  boundary  of  D  between  the  points  of  A  and  C. 

3)  The  point  at  infinity  is  not  treated  by  subtraction  of 
the  pole  at  D.   Instead,  the  problem  is  solved  in  the  inverted 
domain  given  by 

Z  =  1/z  , 

which  maps  the  point  at  oo  into  the  origin. 

4)  The  formulation  of  the  equations  follow  directly  from  the 
theory  so  we  shall  not  bother  to  write  them  down.   We  note  two 
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features.   First,  a  pasting  recipe  is  required,  as  it  is  for 
the  cusped  geometry  model,  to  keep  the  images  on  the  circular 
arc.   Secondly,  the  corner  B  was  treated  in  the  same  way  as 
other  points  on  the  boundary  according  to  the  tangential  shift 
formula. 

5)   Table  VIII  summarizes  the  results  of  the  computations. 
For  these  runs  the  following  quantities  were  used: 

h  -  the  basic  mesh  width  =  l/l2 

r  -  the  fixed  length-to-width  ratio  of  R  =  1.5 

a  -  the  dissipation  factor  for  the  interior  equations  =1.1 

P  -  the  tangential  convergence  factor  on  the  boundary  =  1.0 

a  -  (see  Figure  5-5)  =  .864440 

b  -  (see  Figure  5»5)  =  -4  and  .8 

The  global  difference  between  the  bolutions  for  the  two  values 
of  b  is  that  for  b  =  .8,  the  image  of  the  corner  B  lies  on  the 
vertical  segment  and  for  b  =  .4,  it  lies  on  the  circular  arc. 
The  data  verifies  that  the  convergence  of  the  numerical  solution 
is  O(h^). 
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Appendix  A  -  Tables  and  Figures 

Table  I 
Cusped  Geometry  Model 


Rectangle  length  a  =  Tr/2 

Rectangle  height  b  =  37r/l4 

Coil  radius  r  =  .510075 

Wall  constraint  ^  =  .97 


Basic  mesh  width       h  =  ■Tr/56 

Convergence  parameters  a  =  2.2 

P  =  .9 

7  =  .^5 


Mesh 
width 

Number  of 
equations 

Number  of 
iterations 

Separation 
points  S-.,    Sp 

Midpoint 
M 

a2 

h 

377 

1210 

.24134 

.72657 

1.20279 

h/2 

1425 

2121 

.21851 

.72093 

1.15358 

h/4 

5537 

4187 

.22031 

.71946 

1.14127 

h/8 

21825 

8149 

.21369 

.71888 

1.13635 
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Table  IV 


Contraction  Coefficient  C  (r,h)  -  Plane  Jet 


mesh 
h=l/l2 

r  =  1.5 

r  =  3 

r  =  4.5 

r  =  6 

mesh 
h=l/l2 

r  =  6 

h 

.75506 

.63442 

.62245 

.62150 

V^  h 

.61412 

h/2 

.74891 

.62637 

.61387 

.61271 

V8  h 

.61177 

h/4 

.74817 

.62521 

.61265 

.61146 

V16  h 

.61129 

Table  V 
Contraction  Coefficient  C-|_(r,h)  -  Axially  Symmetric  Jet 


mesh 
h=l/l2 

r  =  1.5 

r  =  3 

r  =  4.5 

mesh 
h=l/l2 

r   -  3.75 

h 

.66872 

.60812 

.60698 

3/2  h 

.60004 

h/2 

.65826 

.5955^ 

.59389 

3/4  h 

.59520 

h/4 

.65662 

.59364 

.59193 

3/8  h 

.59258 

h/8 

.65629 
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Table  VIII 
Conformal  Mapping  Model  -  Convergence  for 
r  -  1.5,  h  =  1/12 


Mesh 

Number  of 
equations 

Number  of 
iterations 

Image  of  the  corner  B 
(Figure  5.^) 

b=.8 

b=.4 

b=.8 

b=.4 

h 

247 

458 

414 

.86444  +  .731371 

.85523  +  .525841 

h/2 

925 

1059 

862 

.86444  +  .722641 

.85565  +  .522971 

h/4 

3577 

2267 

1812 

.86444  +  .720481 

.85575  +  .522271 

h/8 

14065 

4807 

3810 

.86444  +  .719941 

.85578  +  .522091 
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Figure  2.1  Cusped  Geometry  Model 
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Figure  2.2 


Riabouchinsky  Flow 
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Figure    2.4  Vena   Contracta   Flow 
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Cusped  Geometry  Conformol  Mapping 
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Figure  4.6  Treatment  of  the  Separation  Point 
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Appendix  B 

On  the  following  pages  Is  given  the  Fortran  IV  program 
for  the  Riabouchlnsky  model  as  coded  for  the  CDC  660O  computer 
at  New  York  University.   Also  Included  Is  actual  output  for 
an  axlally  symmetric  run  corresponding  to  r  =  1.5^  h  =  K'/l2. 
This  case  Is  set  by  the  parameters 

NH  =  18 

NV  =  12 

on  the  data  card.   Also  Input  on  this  data  card  (Format  14  of 
the  code)  are  the  various  convergence  parameters  and  control 
devices,  each  of  which  can  be  deciphered  by  close  study  of  the 
code  and  the  output. 
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I2?5)  ,PSIV<2?5l,CON(50).       «(    2000)     ,  v  (    2000). PSH     <:Ona) 

UU.  VD 

*L    PARAMETERS 


IF<NW 


1.»D2.MD3.MD4.NC.M.nE.NEP.IP.1S1,1S2. 

,CON(«).FPS 

2,5F9.01 


NE=NE"NV/1?   $NL"(NL 
QstXBl-PI»NV/NW) 
XKP  =  «»<«6L0G(  1  ./ni/P! 
ORT  [  1.  -XKLi'l 
PI/2.  . 1.-KKL21 


Xi(L-'" 

fcpscr 

NVl 

SRI 
CUV 
HLb 
wiNT 


jNTf) 

SCA" 

$R1  =  C 

-XKL2 

IhJT 1 AL  Data 


T=NN1»NV1 


l.-S 
.5/4 


(1)*M 

XK  )/XKL 


tNC«NC*NV/l2 

SXK"f  KH  fj) 
JXKL2=X«i.»2 
SlSl«NI T/151 
JNMiaNH-l 
)S(i";./tl.*C0N(4)«H)SSRIai./Sft 
lACUsi  .'i/k'  $BDVa-?./H 

iHi  =  CnM2)*H  4B3«i;0Nt3)«H 

iUMFsKXL/d.-Xfi.P)     $ERLlM=i.F-h 


I  »  I  TCifi 

■1)*NV/1?*1       5!S?"I52*NV/l2 

S0P»6XP(-P1»NH/NV 
iXKL^FxXKXx.n  ) 
IF«FF (?1/?. ,XKL?) 


call  spe-.'(i) 

GO     10     (27,11 

«J2«1    US£    ^.X 


.98) ,MD2 

rT    PLANP    SOLUTICW    FCB     IMTUl    GUVSS 


i;Q    "•'1    L«i.NHl 

S'JsSKCt  (L-1  .  )«".^.XK.XtLl 

C>J     rS.lWTt  tMAXKO  .  .1  .  "SN     ••i^)) 

(i'^sSOSTd  ,-{XKL*SM)*»4;) 

.>«l!l»(L-l  )»N'V1 

uLI     ^1    NaN^.N2 
1F<.J.F0.NV1  )     GO    TO    21 

S^f  ((N-N1.  )»M,QP,XKP,XKlP1 

5iJPT(AM4yifC..l.-SMP"*2() 

SrjRTi  1  .  -  (X^L  'i»SNP)**2) 

WOLXt (L-1 . )»M. <M-N1 )*H) 

CPi-^UF)  <PF*L(Zi!).fc.XR.O)  ,4lf  AG(Z?1-6F1(AIMAG(Z2).EP.XK^,C*P1  ) 

PLXtXhu^•SN.CN•D^J•SMP••2.SNP•C^F•C^P•D^■••2)*(0..1,)•XKLP•(!JN• 

l)NF-{n.,l,).yKLi:•SN•c^•s^p>l/^r^p••^•xKL2•(SN*S:MP}••^^-(l.-l'KL 

?*C-1P_X(0.  ,EP-*'<L<;-XKP))/XKL 

.'lT.N-:*tL-l)*('^Vl*l)»l)        /li?l-l./Lil'wF/(22-XKP»(0,  .  1.  )  ) 
XIM)B-)Eti   l?l) 
Y  I  M  ■  A  1  U  i  (  Z 1  ) 

iJONTHUf 
X  ( rj  V 1 )  =  C  . 
Y(Ntfl)=D. 


cnp  = 

ONP  = 


XC.NP*: 

X2)* 

|F( 


KD2  = 
1    CAL 


f     IhTr^POLATf    PHtSfcNT    SOLUT10»J    FC»     INITIAL    GUFSS 

.  =»6»^lV(x,Y.zf^,^s.^Hl,NVl,^TCT,xol 


•••  MD2»3  USE  eOUNDARV  DATA  STORED  ON  TAPE  1  tuf  INITIAL  GUESS 

98  REWIND  1 

HEAD  (1)  ZEi.ZS.xn 

00  37  Nai.NTnT 

XlN)BO. 
97  Y(N)"0. 

DO  "**>    L«1.NT0T.NV1 

ZleZf-d'L/NVJ) 

x(L)»REAlIZi) 
96  VtL  )=AHAG(Z1  ) 

N1»MT0T-Ni/ 

DO  -JS  L"N1,NT0T 

Z1BZS(1*L-N1) 

XlL>"«FAL(Zl) 
95  ¥(L)«A1HAG(Z1} 

DO  99  L«1.NV1 

99  X(L)»X'l 

•••  SET  UP  INTERIOR  IMTI*L  GUESS  USING  FIXEC  bOUNDARY  DAT*  FROM 
OPTIONS  FOR  MD2 

16  00  17  L-l.NTOT 

17  eSllL)»0. 

DO  42  L-l.NITP 

42  CALL  ZCHUf;(X,  Y.P<!I  ,KKP.MD4.  (MD^-D'-P.lM'wF) 
iF(MD*.EQ.O)  GO  TO  "3 

no  41  L*-I.N1TP 
41  CALL  ZCrlUGCX.Y.PSl  .XKP.MD4,l,ulNr) 

43  30  5    L=i.NHl 
!Jl=l*(L-l  )*NV1 

5  ZB(lJ»CMPlX<X(N1.).  YIM)  ) 

DO  ."H  L  =  Nl.NTOT 

K=L-M1*1 
3  ZS(t<)sC1PLX(X(L)  .  YlL  M 

iF(MD3.EC3.21   GO  TO  dti 
;•••  PRINT  INl'-AL  HOD'JDABY  VALUES  ANC  tt-BIV*TIveS 

CALL  JtR!  VtUD.VD.PSlV.X.Y.PSl  ,ZB,Zl=.l  ) 

CALL  SPFw(5) 

CALL  iJfcyn/iur,vn.psiv.x,Y,P5i  .Zfl.zs,?) 

CALL  SPe-M) 

CALL  SPE.J(9) 
c?  CALu  SPew(2) 
:•••  BEGIN  MAIN  jTPRATION  iCMEHE 

DO  7  Isl.NlT 
:•*•  THANbFEB  PRESENT  HOUNDaRY  D^Ta  F HO^  ZP.  7«  TO  K,Y 

DO  9  L=l  ,NTOT , NVl 

Zl»    Z8(l«L/NVl) 

X(L)«HEAL(Z1) 
i  Y(L)"AIMAG(Z1 ) 

NlsNTQT-'UV 

DO  9  LbNI.NTDT 

Zle    ZS(1*L-N1) 

X(L)bREAl<Z1) 
9  V(LI«A1HAG(Z1) 
>••  ITERATE  INTERIOR  .  TEXT  EOUaTIONS  4.13*-4,13C 

IFd.EO.IPl.OR.I.Lt.SOO)  |TR«l 

CALL  2CMUG(X,Y.PSI.XKP,MD4.  ITR.UINF) 

IFMTR.EQ.1.AND.I.QE.5001  1P1  =  IP1-IF 

PAGE  NO.  2 


lTft  =  0 

»VG2«0. 

AMAX^O . 

BMAXbO , 

CMAXaO . 
C*.«  ITERATE  FIXED  BOUNDA 

CALL    DERI7(UD.vri.PS 

CALL    EXTR6P(C0N(5) 

00    in    i.s2,NVl 

If  (-10D(  I,  iS2).EQ,0) 
10    Zb<L) =CMPLX(HLENGH, 

2S<U  =  ZS(ll-'B2/Sf3HT 

Ztl{^JHl)«7SCl) 

IF  {■iOU<  1  ,  ISli  .EO.O) 
C«»«     ITERATE    fBEE    KOllNnAR 

CALL  UERiUtUU.VO.PS 

CALL  EXTRJP(CON(M. 

AVGl«(  1  .-EPSl«CO^(e> 
C"*»  CHECK  FOR  lIVFRr.FNCE 

IFUOO.  .LP.RtAL(UD( 

00    13    Lsi.NH 

AvAL»REAi.  (UD(L)/VD( 

HVAL«(REAL(VD(L) )• 
XMjj4)/AVQl-l, 

CALL  MAXVAL( AMAXi *V 

IFMODf  I.IS2).E0.0) 

13    ZR(L l=7H(L)*fi2    •AVA 

C«»»     ITtNATE    L-SHIFT 

DO    i5    L»1.NL 

XV=«EAL( VD(L)«CONJR 

XUaREAL(iJD(L)*CO^JG 
15    AV32sAVG7*(XV-XU)/( 

AV52»AVl,2/NL 

X(jsxt.-hl»AVG? 

lie    ^?'    NbI.NVI 
2H    X(N)»XO 

Ze<l)=CMP|,XlXO,  AtMA 
€•••  CHECK  FOR  CDNVERGFNC 

iF[ AHS( AMix ) .GT.PRL 

CALL  SPEWt?) 

CALL  SPeWl4) 

CALL  Dt»lV(Un,vn,PS 

CALL  SOfc-ilS) 

GU  TO  20 
30  If  (.-^'ODlI.ISD.EQ.O) 
24  iFlMODd,  1S2).E0.G> 

lF('iOt)(I.500).NF.O 

RE-iMJ  1 

-»1TF  Cl>  ZB.ZS.xO 
7  CONTINUE 
20  CALL  SP£wi9) 

CALL  SPEw(lO) 

CALL  SPEwiai 
IF(--iUl.Je.2)  Go  TO 

fleWI^D  1 

-RITE  U)  ZB.ZS.XO 

HO  TO  23 
12  CALL  SPE-'M) 

CALL  SPE-l(5) 

(iO  TO  23 

en:j 


RV  -  TEXT  EClUATICN  4.13D 

IV,X,Y.PSI.Ze#ZS,3) 

XX.NeP.H.NVl,UC.PSIV.C0N(3),C0N(4),l.MD4,ZS) 

Call    MftXUALCCAX.fiEALCVDtLJ/UDtLJJ.L.MAXC) 
AlMAG(Zii(L)l-P5-AlhAG(UD(Ll)) 
2. )«CMPLX(0..AlPAG(yO(l))-AlMAGlUDIl)) » 

CALL  SPEU(3) 
'  -  TEXT  FOUATION  A.jlE 

I V ,  X ,  Y ,  p  s  I ,  <:  e .  z  s .  7 ) 

X«.NE,H.NHl,VD.P5IV,CnN{l).C0N(2),2.MD4.ZH) 
)*EPS»C0N{5) 

I) )»A1MAG( VD{1 ) )  .Oft  .  Y tl J  .LE  .  . 1 )  GO  TO  12 

L)) 

2*AIMAGlVD<Ln»«;)/PSlVIL)»"2«AlMAa(ZB(L))«.l 

AL.L.HAXA) 

CALL  MAXVAL(P^  AX.EUAL-L.'^AXB) 
L*UD(L)*B3  •BVAL*v[:(L) 
xT  EQUATION  4.13F 

( VD(L)  )) 
{UL)(L)  )) 

XV-XU) 


GIZB(l) } ) 

E 

IM)  iJO  TO  30 

1  V  ,  X  ,  Y  ,  ^  S  I  .  Z  B  .  Z  S  ,  I  ) 


CALL  SP6l-<4) 
CALL  SPEM(6) 
OR.  MDl  .NE.21  GC  To  7 


SUflRnUTIvjF     DFRIVluD.VD.PSIV.X.¥,PSl,7e.Z"^.«) 
C*««    COMPUTE    DERIVATIVES    ALO'^G    FIXFD    B^LNCaSY     (K^l    OB    -i)     AND    ALu«li 
C  FREE     aoUNIAPV     {K"2> 

COMMON/m /ADv,PDU,CDV,XO,NV,WH,NVl.NMl,\ToT.vc,S01.CA,Sfll,H,NlT 

COMPLEX    'ID.  vr<,79.7S 

01  MENS  ION    (.iD(ll,VDtl).?8(l).7S(l).x(U.Y(l),'»sI(l  ),PSIV(1) 

GO    TO    (1.2.1).K 

1  30    3    L=1.NV1 
N1»NT0T-NV1*I 

iJD(L)»in"»7'!(Ll*aDV*CM0LX(X(Nl-NUl  ).  V(^'l-NVll  )  ♦CDV«CMPLX  t  X  (  Nl  *?«MV 

11)  .  Y(N1-7.WV1   )  1 

IFCK.En,'')     P<;iv{L)»Anv*PSMN1  )*opv*p(;ifM-Nvi  )*Cnv*PSI(Nl-2*NV1) 

3  vD(L)=':i)"«(Z«;(i  •i)-z*>(L-i)) 
VD(l)  =  <Z<:t7)-Z<:(i)  )-Cnv«2. 
|JD<l)s(Z<;(l)-Zn(NHM*CDv«2. 
\/D(nJ'/l)B'-nu»fC'^NJGCZ'^('^yi-l)l-7S(Nvl-l)) 
RETUWN 

2  no   4    L  =  1.'JH1 
Nal*     tL-^  I. NVl 

•yD('_)s-A"V"ZO(i   J-qliv-CPl.XCXfN'U,  Y(V4i  )  1  .COV»rMR|.X  (  X  (  N»2  )  ,r(N*2)) 
3S|VtL)s-APU*P«;I(Nl-»D^»PSl(^*l)-C'^v«F<;i(N*2) 

4  'JD(L)=CDV«(ZO(i  *t  )-Z'^  (L-1  )  ) 
JD(l)=':nu*(ZR{3)-C'^PLx(2.»X0-BFA|(7Pf2n.AlM4G(Za(?l)l) 
'JD(NHl)atZB(^lHO-Z^('^H))»C'?V*2. 
VDtNMlJstZSO-ZSdjl^CCV*?. 

RETURN 
END 

SU'^RIUTpiF    fi   1D(X,  Y,R,'*,RFS.  X'<P,M,K  C  ,K».N  V  ,L   I^F.HP4,  ITH) 
C»»-    SFT   aouNHARv    rOwnITI0^S    ALOvG    I^Te«;Fif:E      "s-l    R&GtO^'   »    TO    REGION   oi 
r  Ka-l    RCG!0»l    Rl     TO    RERION    R 

nl-E'jsiov  x(i  ).y(i),P(i) 

'.-O-^PLEx    7 

N25Nr-i 

DO    1    L>1.NA 

IF(L.EO.MA)     w2rNl 

DO    2    ■yeN'  .N2 

IFdTR.Cj.l.ANn.K.Ef!.-!)     P(%)sP{K)-LlNF/(l.*''n4).vtN)**(l»MD4) 

r.K     TMOLXf  X(Nl,Y('«))*'<«BES/rMPLX(  f  L  -  1  .)••■.( '-C-2  .  •'  ♦»J-^l )  •H-XKO  ) 

YUjoAlMaltZ) 

IFUTH.FT.I  .A'jn.K.EC.     1)     P  <  K  )  sP  (  v  )  *(.  I  ^  F  /  1 1 .  ♦"O*  ) -Y  r  N>  ••  ( 1  *HD4  > 
?     jCIV)sUFAi   (7) 

VIlaNl  -MVI  •! 
1     N2«*J1.»1 

uETiJRN 


"lUflROUTI^'E     MAXVAL(XMAX.XW«i..'«.MAV| 
C*»«    CALCULATION    OF    *HXIHuM    MOUNDARY    RFSIfALS 
IF(AaS(X>'AL  )     .I.T.     *B5(XMAX))       RETL»N 
XHAXbXVAL 
MAXbN 
ilETURN 
^NO  PAGE  NO.    t 
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SU9R0UTINR    ZCH|]G(KtY.P.XKP.MO«.  ITB.U  t^F) 
(;•••     ITtRATE     UTEfllOR  MD«  =  0     PL^NE  ITR«0     DO    NDT     ITEfiATfe    PSI 

C  MD*  =  1    *XUL    SYM^-ETPY  ITRil     ITERATE    PSi 

(;0«rtO^/Sl«'*DV,BT)V,COV.)(O.NV.NH,tjVl.K''l,MOT,NC.SBl.C*,SHI.H,NIT 

OMcNSIO^   x(l),V{i),P(H.Yi(?nO0> 

NAB-gVt-nC    I    NA1  =  'J**1    S    N»2  =  N4-1     $    Nl»<    5     NC1»NC*1     $    RES«l./OlNr 
C»«*     ITtHATe    THBOUfiH    Rt-GION    R 

DO     1.    L53.NC 

1F(  ITR.EO.O)     GO    TO    1 

IF  t  ^')4.i:.7.0)    P(|.>=Sfll*P(L)»t;»«(P(L*l)*P(L-l)*2.*P(L*NVl>) 

irt  .HD^.C-J.l)     P(L)sSfil*PlL)»(P<L*l)/<V{L>*T(L*l>)-PlL-l)/("l'(L)*Y(L- 
XI)  )»2.»P(!.-NVll/(Y(L'*''(L*hVl  1)I/(SPI/(V(L)*Y(L*1)  )  "SR  I /(  Y  (  L  )  ♦  Y  (  L- 
X1JW2.*SRI/(Y1L)*Y(L*NV1))  ) 
1    Y(L)aS»l«''(L)*r4-(Y(L-l)*Y[L-l)*?.*Y(L*KVl)  ) 

JO    11     LBp.NH 

IF(  ITfl.E'i.O)     (iO    TO    6 
IF(HD4.E0.1)     GO    TO    7 

no  6  N=-ji-'y2 

6  P(\)*SRl«PtN)*rA*tP(N*l)*P<N-t>*P(N*'iV]J*P{N-NVl)) 

t»o    ro   6 

7  JO    «    >,="J1.H2 

9  P(\)BSai«P(N)*tPtN*l>/fYtN»»Y(N*i))*PlN-l)/(Y(N)*Y(N-l»)«P(N*NVU/ 
X(Y(N)*Y(N.Ntfl)  )•P(N-^JVl)y(  Y(ly)*v(\-hVl)  )  )  /  (  SP  I  /  (  Y  (  N  )  •Y  ( f4*l  >  J»S»I/I 
XY(N)»Y(N-in*SRI/(Y(h)*YtN*NVl>)*ERI/tY(*J)*Y(N-fjVin) 

■S    JO    11    \s\l,HZ 

<(N^«SRl•x(^l)*CA•lX(N•l)•x(^-l)•x^N•^v1|•x(N-NVl)) 

11     YtNlsSRl"''(N)*CA«lT(N*ll*Y(N-lU¥(N«VVll*T(N-MVll) 
IFt^.Lt.NA)     GO    TO    IC 
4SS2-1 
X(>\)B5.*l-XtN)*CA"(2.*X('*-l)*X(N.NVl)*ll\-NUl)  1 

-•••    C>*OSS*J«e:?    to    PgRION    HI 

':aH.    '^Llo(X,Y,P,-l,fifcS.X'<P.«,NC,N«J  ,^Vl,'I|f»F,MD4,ITHJ 
U  (  'TR.Ef).l.ANn,MD4.fcQ.l)     CALL    CRCSSt     l ,  X.  Y,  Vl.P,  HES.H.XKP.  vv  ,  Nn. 
X.JA?.SC.'^''1  1 
C««»     lTE*taTt     THHOur,H    RPGlON    Kl 
yO    1?    LaNCl,HV 
If  (  iTfl.ET.O)    'lO    TO    li 

|F(Mr!4.S3.CI)     P(L)=SHl«P{L)»C*»(P(L*l  )*P(L-1  ) -Z  .  •»  (  L-NVH  I 
IF(1D4.eo.  I)     P(L)=S«1«P<L)*1P<L*1  )/(Yl(L)*YltL»in*PtL-l>''(lfltl->*^ 
XllL-t))*2.«P(L-NVlI/(Yl(L)*vi(L»NVl)))/(SfiI/(YlCL)*Yl(L*ll)*SRl/(Y 
Xl(LJ*ri(L-l J)♦?.•SRI/(Y1{L)•Y1(L*^V1  )  )  } 
13     Y(L)»5»1«''(L)*CA»(  Y{L•■l)•YCL-l)•2.•Y(L•^Vi)) 

N1  =  N.-1 

UO    1*    L«?.NA2 
n=Vl*NVl*l 
V2  =  ,\?»MV1 

IF(  iTi(.EO.O)    GO    TO    1> 
lFCia4.eo.l)     GO    TO    16 
tiO    l7    ^J•V1,N2 
17    HtM)sSt*l*P(N)*rA*(P('»*l)*P<N-l>*P(^»'^Vi)*P(N-NVl)) 
.iO     TO    15 

PAGE  NO.    5 


16    00     !■*    NeM,»J? 

IS    9(M)»SRl«P(M*rP(N*l)/(Yl(M)»yl(...i  )  J-Pfg-I  )  /  t  Yl  (  N  )  .Yj  (  N-i  )  )  ^p  (  n»N 
XVll/tri(M.V1  CJ^^VDUBt^J-NVI  )/(Yl(K)«¥i  tN-NVl))l/{SRl/{YHN)*vi(N 
X*l)I*SRI/(Yl(»4)*Yi(N-i)).sqi/(Yif\)*vi(K»NVin«SRI/{Yl(Nl*Yl(M-NVl 
X)  )  } 
IS    DO    Ifl    «J»M1,W9 

K(^)•SRl•X(Nl♦r*•(X(^•l)•Xf^-l)•lr(^•^JVl  )*xtN-NVi> ) 
19    YiNI'SRl^ftNl^rA-CYCJ.D.VfS-D.YC^^wVll-YCN-NVm 

N«N?*1 
14    X(N)«SRl«X(N1*r*«(2.»X(N-l)»l((ty*^vn»X(K-NVH) 
»JbNA»NV1 

XCJ     i=SRl*»l'I     )*Ca*(?.«X(M    -1)»X(K     •NV1)*X<N    *NV11J 

1F(  ITR.E0.1.4»in.Mo4.eo.l)     CALL    CPCSS(-1  .  X,  V»  Yi  ,P  .  RFS  .  H.  XXP  .  NV  ,  »VA1 . 
XM*?,MC.NV1] 
C«««    C^OSSnvER    TO    REGION    R 

CALL    fLlo(X,V,P,     l^RfS.XKP.M.NC.*!*!  .^iVl  ,  JI'WF,hD4,  ITRl 
HETUMN 


SUPHl:lUTI^E    RP*nYtx.Y,Zti,2S.VHl.\Vl,\TCT1  .l"! 
r»»«     INTFRPOLATE    SntUTIO**    ALONft    POPNDAPY    FOG    Pi^*'    Q^F  If-'^t^^HJ 
r  FACTOR    OF    2 

TIMP'^Sin'y    )fCi  l,Y(11.7R(l),?S'l) 

COHPLFX  7S.ZP 

DO  515  "Jal  .NTOT1 

X  ( 'J )  s  C  . 

50  vcg).c. 

;iO  51  l.*t  .NH1  ,9 

Nl  =  l*(l  -1  )*NV1 

X{'Jl)=^i;4L(ZR{l*L/2)) 

51  YtMDsAl-AGfTSd'L/Zl) 
10  ^?    Lai  .NV1 ,9 
Nl*>JTaTi.L-^Wl 

V(  Jl)=RE»l  fZ';{1*L/2l) 
5?  VCM1  =  »JN1R{7S(1*L''?1  1 

KCI  =  X(11 

00  53  L*1 .NV1 
5:^  X(LI«XO 

■)0  4fl  Ls5,»JHi  ,? 

Nl  =  vjVl-fi  -?)*l 

'{Wlle.l?5*(.X(Nl-J»^'Vn*6.»t(Kl-NVJl«?.«X(M*hVM) 
4n  Y(^'•,  )  =  .175•r-Y^Nl-i•^Vl>•6.•V(^l.^vI  )«"»,*YIN1«MV1  )) 

f)0  41  I  •«:,»; V1*9 

Ml»^'T^Tl*L-3-NWl 

X(Nl)  =  .135-(t.»X(Nl-l  )*A.*y(NI*ll-X(''l*7)) 
4]  Y(»Ji)s.i?5«(T,*Y(Nl-!J*6.*Y(N.i»U-V(*.l*!)) 

^is^*-!*! 

X(V'_)s.S«{X(  N1*NV1  )*X(  M-\V1)) 

Y<11)s.5«(v(wi,NVi  WY(  Ni-^V1)  > 

IJiVTOTl-l 

it(M)sX{VTOTl) 

YtNl)3.5**tNT0Tl-?I 

.RETURN 
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SUURDJTIVE  SPEWfNSPEW) 
C***  DIVULGE  ArjSMERS 

C   NSPEMbI  initial  data  «  ITERATION  DATA 

C         2  HEAOfcR  FOR  ITERATION  DATA    7  CO\VERGENCE  DAT* 
C  3  FIXED  BOUNDARY  DaTA  6  \OCaL  LINES 

C         A  FREE  BOUNDARY  D*TA  5  FSI  CERIV  AlOnO  FRgE  BOUNUaR' 

C         5  DIVERGENCE  OaTA  10  HRAr,  CALCULATION 

COMM0M/8l/ADV,RDV,CDV,X0.NV,Nh.NVl.\wl,KTOT.NCiSRl.CA,SRl.N.NlT 
C0MM0N/8?/Xr(L.XKLP.XK.XKP.E.EP.0.CF.  ISl,  IS2 
l.<lDl>-HD2.MD3.HD4,AVGl.AVG2.AMAX,eKix,CKAV.NAXA.HAX8.UAXC,Bl.t]2<33. 
2SR. i .ERLIM'HLENQM.UINF, 1  P. NE. NEP, Nl . PFS. ?B t 225 ) . ?S ( 100 ) • 
3UD(225),vn(225),PSIVl225).CONl50).   » (  2000)  .Y(  20uO).PSI(  2000) 

oiMt^NSiON  NSdnoJ 

COMPLEX  ZB.2S.Z1,Z2.EZ.UD.VD 

GO  TO  (1.2.3,4,5.6,7,8.9.2%),NSPFh 

1  PRHT  10.XKL.SR.CON(4),M,B1  .  CGN{1 )  .  X«LP.  Nl  T.  NMl  .P2.  C0N(2)  .  XK.  ERL  I" 
l.NVl,33,C0N(3).XKP.  I  P  .  NTOT  ,  NE  .  E  ,  MC  1  ,  MLEKQH  ,  NEP  .  EP  .  MD2,  U  I  NF  ,  t-'S  .  U 
2.MD3,NC»0P.HD4.NL 

10  ►"OrtHATClril.  9X.•■<L•.^1^.8.4X,•RELAXATI^^  F  aCTO**"  .  2F7  .  4  .  7X. 

l»MES-<  S I  Z£  •,  ri  1.6, 5  X."  CONVERGENCE   B3*,5F7.4/5X..«LFRIME*-Fl2.6.3X 
2.-MAKIMJM  ITCRiTIONS".  114.  7X.»H0fiIZ  FTS,  .  J  ir  ,  9  X  ,  .F  ACTORS  Bf. 
3^F7.4/H)t,»K»,Fi2,8.1UX,«ERR0R  L  t  ^' n».  f  14  .  ?,  8X.«VE»T  PTS».I1<' 
4.18x,-d3».2F7.4/6X.»KPRiME»,Fl2.a,]«»..!o», I1«.1X, 

5«T0-iL  EQUATIONS*. ll0.4X.*PTS  1 N* . ex , •WE* .  I  7 /ll X , •£• . Fl2 . 8. B> . • I NP 
6UT  SET  MCI1«,  114. 4X. "CAVITY  lENGH*  ,  F  1  ?  .  ft  ,  4  X  »  •£  XTRaP    NEPRIM-*,!?/ 
76X.«EPRlMe.,Fl2.8.18X,»MD2». U^.IX.-VEl  »T  iNFlNITY^.riO.b.JK.^cXI 
8RAP  4r  6PS«.F7.2/llX,»Q#,Fl2.J«,!ex.«*'03«,  1 14,  30X  .  •  |NTE«F'>-5 

9  NC». I7/6X.*0PRIHE«.fl2.e.l8X.«Mn4.,  1 14 , 32X , -L-Sw 1  FT • , 5X . -NL • , 1 7) 
RETURN 

2  PRI JT  13 

13  F0RM4T(///5X, •ITERATION  DAT  A*//2*.  •  I  TER      AMAX-INDEX      imx-IN 
1D6X      C^'AX-INOEX      MDPT    SEPF  T     1/OOSO     1/OlSO     1/02SQ 

2  SLJPEl   SLOPE?     RES-XO         XD»/) 
RETuflN 

6  Zl»CMPLX(*|MAG(ZRt       I »  )  # A  I MAO( ZB ( K"l ) ) J 

PRINT    14,I,AMAx,  ■4AXAiaHAX,HAVe>CMAX,MAXCtZli  AUGl.CONd  },CON(  j) 
1    .C0N(2).(;0N{4)  ,»va2.xo 

14  F0RMAT(l)(,|S.3{E11.2.  l4),F9.9,Fe.5.2F10.*.2F8.4.2FiO.S) 
RETURN 

3  PRINT    25, I . (ZS(L).L«l.Nvl) 

25  FOP  itT(//5X.*f"ixC0  BOUNOaRv  aT  1  TEPaT  !  ON*.  I5//(  7  <F10  .  6.f  9  .  6  )  > ) 
PRINT  19.  ((un(L).VD(L)  ),L  =  1.'»V1) 

15  F0RMAT(//5X, •FIXED  BOUNDARY  DERIVATIVES   XU/YU/X V/ Y w. , // ( 4 ( ^X . 4t 6 . 
X5)  }  ) 

RETUR.^ 

4  PRrjT  la.  !  .  <  ZB(L1  .L«1.NH1  ) 

IB  F0R1AT(//5X..FREE  BOUNDARY  AT  I TgG A T t CN« , I  5// ( 7 ( FlO . « , F9 . 6 ) ) ) 

PRINT  17,  t (UD(L).VDIL) ).L=l,NMll 
17  F0RHAT(//5X*»FRE£  BOUNDARY  DERIVATIVES   XU/YU/XV/VV* , // ( 4 t 2 K . 4C fl , 5 
X))) 
RETURN 
3  PRINT  19.  I 
19  FDRMATI////10X. .DIVERGENCE  AT  1  TER*T  ICN«  ,  15  ) 
RETURN 

7  PRINT  20.ERLIM 


2u  FOHMaT(//10X..PRR0R  is  BELOM^.ElO. 1.// ) 
QO  TO  6 

8  PRINT  85 
N1»2-NV1 
n2iHC*1*'IV1 
DO  81  L32»NH 
00  (i<f    Nbn1,n2 
\3  =  N-'J1*1 

XU«(X(N»1)-X(N-1) >«*2*(YtN*l) -Y(N-l))^»2*(X(N*NVl)-X{N-NVl))"^i' 
l-tY(N*NVl)-T(N-NVlJ)*«2 

tfCXU.GT.QJ  NStN3>»l 

IFUU.LE.O)  NS(N1)«P 
8^  COMTINUc 

PRl  JT  83. ( NS(K) ,<il,N3) 

NlsNl*-iMVl 
81  N2S1IN0  t'-lC-L^CNVl^D.L^NVl^NV) 
87  FORHAT[//iox,76HLOCaTION  OF  NODAL  lInES  OF  ( x  V»«?*Yi/«»2-XU*»<?-YU»^ 

X2)    1   IF  (*)      0  IF  {-)    //) 
1?  ^0»^AT(5x.lOOI1 ) 

HETuRN 

9  PRINT  16, tPSl VtL) .L^l.NMl} 

16  fowiAT(////irx,»»siv   ALONG   FflEE   BCu^OASY•// 1 IX.  i<;fo.  *)  ) 
R6Tu-*-g 

2i    CALL    ueRIV(un.VD.PS]V.X.Y,PSI.ZB.ZS.3) 

PRI'fT    27,  {PSIVtL).L»l.NVl) 
?7    FOSHA"(//10X.«PSIU    ALONG    FIXEC    BOL^Z AC Y«// ( t X. 15F9 . 6 1 ) 

CON(l )=J-ftL( 7a(NHl) ) 

C0N(2>=AIM4GtZ«{ll ) 

C0N(i)»AMA3{ZR('JHl  )  ) 

C0N(4)sSTRT(l ./AVGl) 

SGtX«i.»XXLP/(l.-XKLP> 

SGCmP=1./AVG1/UINF-«2-1. 

CD=XB/>.^(EP-XKL««2"X'<P)/(XKLP»«2*t'^-XKL«*2*XfP> 

ASO. 

d  =  0. 

JJO    21    L«3,NV,2 
21     AsA.oSIV(L)^«2/AlHAG(VDtL)>/AIMAG{ZS(L)>««HD4 

jO    2?    L«2tNV.2 
4J>    Hsq-PSI  V(L>^«2/AlMAG(VD(L)  1  /  A  I  M4G(  ZS  (  L  )  >^«H04 

CDC-MP  =  l  .♦(l.*Mr)4)«AVGl/AIMAGtZ»(M-l))«»(i*HD4)«rM/3.*(4,«B*2.^4)) 
PRf-iT    2fj,_0N(l  J,C0N(2),C0N<3)  ,UI\F,::iNi4).SG?X,br,c;MP,C0EX,C'JC*fP 
.»fi    FOR-^4T(//lOX,«nR*G    DATA«//l6X.»SEFAKeTCN    POINT       »•  .  FiO  .6.  1 'X. -CAV 
XlTY     HcIS«T..F10 .A/34X,»Y*,Fi0.6,l3x,.SPEPD    AT     I \[ ' M TY« » F 10 . 6/I6X . 
X^F-lEf    ?lOiJNDARY    SPEED^.fl0.''//5X,«EX4C'     V4LUES«»1JX,«SIGMA*-'^10.6.3 
XX,«CO.-<PUTf-D    V4LUPS^.7x*«SIGMA«.Fl0.e/3»  .•(PLANE    CAS€    ONLY)  CD/(1 

X-SlG''4)*.F10.6,l8X,«CLVtl*SIGMA)^,FlP.iS) 
RETuR-J 
END 
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SUfROUTIVE    COCStK.K.  V.  vl.P,B6S.U.KKP.N«.N»l-,N»2.Nr.NVl) 
-•    S?T    V    ViLMES     IN    RfeGlDM    **!     4ND    PSKSTiP)     ALONQ     lNTPRr»CE 

OI-^ENSION    Xtn.VdJ.PtU.YKl  ) 

COMPUEK     71, Z> 

NISNC 

DO    1    U'l.NAl 

N2=N1*1 

IFtL.EO.M*l>     N9"Ni 

00    2    NsN' .N2 

Zl»l.^(CMPLXn((N)-X(l)     ,V(N))*fl6«t/CMPLXt(L"ll«W,f'yC-2,*L*N-Ni).M-K 
1<P)) 

IFtK.eO.')     Yl<N)"*HaQ(71) 
?   P(N)«Pt'gi«S0PT(PeiL{2l»C0NJG(Zlin»«'< 
'.     NleNl*NV»2 

|F(K.EO.-U     OETUPN 

M1"NC*2 

■J2«*Jtfl 

DO    3      L=1 ,M49 

00    ^       N«M] ,N? 

IC"(N.60.'"V1I     fin    TO    4 

Vl<'y)a4lMAR(l  ./(CMPLy()((V»-«(l)      .  >(  K  )  )  *R  "=  C/CPL  X  (  (  L  "1  )  •  H  .  (  N  -  n2  .  N  V  ) 
1»M-X<P)  J  ) 
4     CONTINUP 

3    N2«N2*NV1 
VKNVl  JO. 

Rgf JRN 
END 


SuRRnuTIvc  EvToAPfrq.x-.'i.u.sjvi.lir.P'^IV.il.Al.^.HD*.?) 

•  E«TfiAPOL»T[nN  rnR  i/oJS'J   (H=2)    a^p  wssc    t^si) 

niMENSIO^    UD(H|PSIVtl),!(\) 

C"nFT(Zl,7?,2^,74)  =  'i.z«-72*Z'» 
AS".     S     R=0,     <     rcD.     ^    D=P. 
no    1    Lsi.N 
(10    TO    C*,?)," 
?    i<»hVl-L 
MA-L-H 

ns0*3EALlJft(lf)*C0NJ3(U'^'(«))WP'=Ivif)»«5«4|M4r.  (;(K))«»(?"MD<1«(-M«L 
1     ) 

50  rn  4 

3  "(eL*! 
A»A«L«H 
D«0*REALt'jri(«)«C0fgJCtuO(K)J)/PSIv(t.l»«5.4I>'Af;(?l««l)«*[?"«J<)»tH*L) 

4  C"C»REAL(UO(Kj#CO^JG(UDlK)  )  1/PSIV(K  M^Ptil^A-^eZC*)  )«*(?«m04) 
1    a»B*(L«Hl"»2 

DET5rOPTfrLOiTfN),4,A,«) 
AO'FOETtr, A,n,n)/DcT 
41»nET{rLO»T(«j),C.».0)/nET 

rs'AO 

RETURN 
END 


^UNCTIflN    PrKu.e.XX.O) 
•    rORMOLA    fllS.OI     PAGE    Sntl    RPrERFNCi:     (4) 
P|»3.141  =  92«53'!897V 

4sPI»U/XV   $^U"=0.   Sp=1.   S0?=1.   JIEeC**? 
DO  1  L*1.100 
Psp.O 
P2  =  P?«':2 
■*»P/(1  .•i'2) 

IF(-II.L=.    l.E-15)    GO    TO    ; 
L    SUMbSiJH*(!.5Ipl.(|  *A) 
y    PFl  =  (2.*o!-S'iM    ♦E.Jjyx'* 
RETtJtJ'j 


ruixcTirTj  crtnp^i.nK^O) 

•    fOR'-ULA     9nj.3f     paGf    2«9    Pt^PfiiF  NCC  f  4  ) 

50    1    L=1.4S0 

Cf'd  .^-'  )  /L*f:F 

"Ts(  (2.,,  -1  .  ,»nT-    >I''t''P<;i  >.»(5.|  .,)♦    rc*t((  PCI  j)/(?.*Ll 

^fcTu-i^ 


fU>':Tin--l    SNFfH. ■).)(<, "tCD 

^rtKWjLA  one. 01   oabe  3113  peFCfic\TP   /4) 

PI=3.l4iK9?6^3q8''7« 
•1=1*34/4,  3B(i  ,/o) 

'^Q    I    1  =  1.  ^J 

SNFsS*jr*f)«.[»SIV(  t?»I*J)»A'/Il.-''.«t?«T.l)) 

=:Bjr«e;NF»<:ooTrOM^.«pi/(XK»¥-L ) 

^ETJR-j 
F-Ml 


•  ••     rnt/'^JI.  4     310, C*!     P»5i=     ?97     ■*CFF*it\i-C      (4) 
*iU"=1  . 
!10    1    L=i.?n 
1     SL"lr«U-*<.«0«*l  /'l  .*'>«-(?»L)l 
r-l53.14tS'S2*S-^St9''»*iu-./P. 
RETji-N 


JUNCTION    rxKi  (yK,Q) 
•  •    P^AMOLE    -^    J46f'479    -*eFF»F';CF     (]H 


rjO 


=1.110 


•('■ 


("*-lU 
l.E-13)     (50 


'0    3 


PAGE  NO.    9 


fftSI     .L= 

fl'''L"3.1«1^996«!3'38»79/!<K*2.»^CftTfCJ«*ll.'-«*2 

^lETja-J 
END 

PAGE  NO.    10 
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This  report  waa  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Commission,  nor  any 
person  acting  on  behalf  of  the  Commission: 

A.  Makes  any  warranty  or  representation, 
express  or  Implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  Information  contained  In  this  report, 
or  that  the  use  of  any  Information, 
apparatus,  method,  or  process  disclosed 
In  this  report  may  not  Infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  Information,  apparatus, 
method,  or  process  disclosed  In  this 
report. 

As  used  In  the  above,  "person  acting  on  behalf 
of  the  Commission"  includes  any  employee  or 
contractor  of  the  Commission,  or  employee  of 
such  contractor,  to  the  extent  that  such  em- 
ployee or  contractor  of  the  Commission,  or 
employee  of  such  contractor  prepares,  dis- 
seminates, or  provides  access  to,  any  infor- 
mation pursuant  to  his  employment  or  contract 
with  the  Commission,  or  his  employment  with 
such  contractor. 
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